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FOREWORD 


The  Multisatellite  Filter/Smoother  system  of  computer  programs  was 
developed  by  the  Space  and  Surface  Systems  Division  of  NSWC  in  support  of  the 
Defense  Mapping  Agency  various  geodetic  applications  of  GPS  and  the  Navy’s 
Stratenc  Systems  Program  Office/Applied  Pnysics  Laboratory  SATRACK  project. 
The  pnncipal  software  developers  were  D.  Clark,  K.  Davis,  E.  Durling,  M.  Eward, 
and  H.  Ball  as  members  of  the  Physical  Sciences  Software  Branch.  This  software  is 
part  of  the  new  OMNIS  ot  bit  computation  system  under  development.  I  am  grateful 
to  MCSI  and,  especially,  Susan  Bowen  for  helping  in  the  production  of  this  report. 
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INTRODUCTION 


The  Global  Positioning  System  (GPS)  is  a  passive,  all-weather,  worldwide  navi¬ 
gation  system  that  utilizes  ultrastable  atomic  frequency  standards  to  provide  navi¬ 
gation  messages  and  signals  of  the  required  accuracies.  The  Space  and  Surface 
Systems  Division  of  NSWC  has  been  tasked  by  the  Defense  Mapping  Agencv  and  the 
Navy's  Strategic  Systems  Program  Office  to  develop  estimation  techniques  for 
highly  accurate  epnemeris  and  clock  determination  (aftei  -the-fact)  for  geodetic  and 
the  SATRACK  applications,  respectively.  For  geodetic  ^plications,  such  as  point 
positioning  and  satellite-to-satellite  tracking  involving  GpS,  the  user  is  interested 
in  accurate  ephemeris  and  clock  information  at  all  times  and  with  minimal  disconti¬ 
nuities.  For  tne  SATRACK  application  (missile  tracking  primarily  during  powered 
flight),  the  user  (Applied  Physics  Laboratory/Johns  Hoj^ins  University)  is  inter¬ 
ested  in  tliis  information  over  a  particular  geographical  area  for  a  short  time  span. 
The  goal  of  this  development  effort  was  to  produce  common  software  that  has  tne 
flexibility  to  optimize  the  orbit  and  clock  estimates  as  required  and  generates  the 
required  products  for  these  two  primary  applications.  In  addition,  tnis  software  was 
to  nave  extensive  diagnostics  and  be  easily  reconfigured  and/or  modified  to  be  used 
as  a  research  and  development  tool  for  accuracy  evaluations  and  other  studies. 

The  result  of  this  development  effort  is  called  the  GPS  Multisatellite  Filter/ 
Smoother  (MSF/S)  system  of  programs.  Smoothed  range,  correlated  range  differ¬ 
ence,  and  two  interferometric-type  derived  observations  based  on  simultaneous 
range  observations  can  be  processed.  A  Kalman  filter  followed  by  a  smoother  was 
chosen  as  the  estimation  technique  for  several  reasons.  The  unmodeled  acceler¬ 
ations  acting  on  the  satellites  (due  to  modeling  deflciences  in  the  gravity  Held, 
radiation  pressure,  and  thermal  radiation  models,  as  well  as  control  system  induced 
effects)  and  ^e  random  behavior  of  atomic  clocks  are  best  handled  by  stochastic 
estimation  techniques,  i.e.,  Kalman  filtering.  Fixed-interval  smoothing  can  be 
accomplished  because  processing  is  after-the-fact;  thus,  estimates  at  a ^ven  time 
can  be  based  on  botb  past  and  future  data.  A  square  root  information  futer/smoother 
(SRIF/SRIS)  formulation  utilizing  matrix  triangularization  techniques  was  selected 
for  Uiis  system  primarily  because  of  its  numerical  accuracy  and  stability.  In  addi¬ 
tion,  ^is  smooming  procedure  requires  that  only  inverses  of  upper  triangular 
matrices  be  computed,  as  opposed  to  inverses  of  full  matrices-as  is  the  case  for  most 
covariance-related  smoother  implementations.  Also,  if  smoothed  covariance  infor¬ 
mation  is  not  required,  the  upper  triangular  matrix  to  be  inverted  is  only  for  a  subset 
of  the  parameteis  in  Uie  selected  implementation.  These  ideas  will  be  expanded 
upon  in  subsequent  sections. 

The  multisatellite  capability  was  adopted  because  it  is  the  optimum  way  to 
separate  the  satellite  clock  offsets  (of  interest  to  users)  and  the  station  clock  offsets 
(of  no  interest  to  users).  This  is  because  of  the  simultaneous  tracking  of  a  giveu  sat¬ 
ellite  by  two  or  more  stations  in  conjunction  with  simultaneous  tracking  oif  several 
satellites  by  each  station.  Intersatcllite  orbit  and  cluck  covariances  are  required  for 
the  SATRACK  project  and  are  only  available  with  simultaneous  pi*ocessing.  Also, 


1 


NSWCTR  87-187 


this  allows  the  system  of  programs  to  accommodate  doubly-differenced  data  types  in¬ 
volving  two  satellites  and  to  eventually  he  expanded  to  accommodate  the  proposed 
cross-link  ranging  (satellile-to-satellite  tracking)  data.  In  addition,  the  adoption  of  a 
multisatellite  capability  affected  the  formulation  of  the  Kalman  filter  process  noise 
incorporation  dealing  with  minimizing  array  storage  requirements,  to  maximize  the 
number  of  satellites  that  can  be  processed  simultaneously. 

The  purpose  of  this  report  is  to  provide  the  assumptions  and  mathematical 
details  for  the  GPS  MSF/S  system  of  pro^ams.  First,  tne  overall  GPS  data  flow  is 
given  to  introduce  the  reader  to  the  preliminary  computations  required  before  exe¬ 
cuting  the  MSF/S.  Then,  the  estimation  concepts  emploved  are  discussed;  followed 
by  the  timeline  definitions;  and  a  detailed  description  of  the  adopted  state  emiations, 
underlying  models,  and  resulting  process  noise  covariance  matrices  rcauired  by  the 
Filter.  Next,  the  observation  equations  and  the  corre^nding  partial  aerivatives  of 
the  data  with  respect  to  the  parameter  set  are  given,  l^en,  the  overall  Filter/ 
Smoother  processii^  flow  is  provided  to  establish  the  groundwork  for  the  detailed 
descriptions  of  the  filtering  and  smoothing  algorithms,  the  generation  of  the  solution 
and  diagnostics,  and  the  propagation  of  trmectories.  Derivations  and  other  relevant 
technical  details  are  supplied  in  the  appendices  to  enhance  understanding  the  main 
text. 


GPS  DATA  FI.OW 


The  possible  observations  are  pseudorange  and  range  difference  derived  from 
integratea  Doppler  or  phase  measurements.  Preprocessing  of  observations  may  be 
rec^uired  for  various  reasons.  Smoothing  of  pseudorange  measurements,  accummu- 
lauon  and  differencing  of  integrated  Doppler  over  longer  intervals,  elimination  of 
duplicate  data,  and  assignment  of  observation  standard  deviations  and  pass  num¬ 
bers  may  be  required.  Tne  pass-oriented  data  must  be  merged  into  a  time-ordered 
format,  and  time-ordered  data  in  one  format  must  be  converted  into  another  format. 
Data  from  various  sources  must  be  merged  for  the  time  span  of  interest.  The  MSF/S 
formate'"^  tains  the  following  information  for  each  observation:  time  of  observation, 
observation  value,  standard  aeviation(s),  editing  flag,  data  type,  integration  interval 
for  range  difference  data,  station  number,  satellite  number,  channel  or  tracker  num¬ 
ber,  pass  number,  source  of  weather  data,  as  well  as  temperature,  pressure,  and 
relative  humidity  in  order  to  make  the  tropospheric  refraction  correction.  Figure  1 
gives  a  simplified  flow  of  the  GPS  data  after  the  observations  have  been  pre 
processed. 

Reference  trajectories  covering  the  span  of  interest  are  required  for  both  the 
Corrector/Editor  and  MSF/S  systems  of  programs.  These  tre^ectories  are  integrated 
either  in  GPS  or  UTC  time  using  initial  conditions  from  a  previous  fit  cr  obtained 
from  the  GPS  Operational  Control  System.  For  the  Corrector/Editor  system  in  its 
normal  mode  of  operation,  the  only  information  required  from  the  trajectories  is 
inertial  position  and  velocity  (obtained  by  numerical  differentiation)  as  a  function  of 
time  and  polar  motion.  However,  for  the  MSF/S  system  the  reference  trajectories 
and  the  corresponding  Lagrangian  interpolation  procedures  provide  the  following 
additional  information: 
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1 .  Homogeneous  variational  equav*ion  solution  for  the  epoch  orbital  elements, 
i.e.,  partials  of  current  time  position  and  velocity  with  respect  to  epoch  or¬ 
bital  elements  needed  to  relate  changes  at  epoch  to  all  other  times.  (The 
cavity  field  model  on  which  the  filtering  is  oased  is  thus  determined  by 
tile  integrator  along  with  the  inertial  reference  frame.) 

2.  Initial  values  for  the  radiation  pi'essure  model  parameters  and  the  body- 
axis-to-inertial  rotation  matrix  and  body-axis  radiation  pressure  accelera¬ 
tions  at  the  required  times.  (These  quantities  are  required  for  the  compu¬ 
tation  of  certain  partial  derivatives  associated  with  solving  for  sto^astic 
radiation  pressure  parameters  and  permit  evaluation  of  the  complex  radia¬ 
tion  pressure  force  model  only  during  the  integration  procedure.) 

3.  Partial  derivatives  required  to  solve  for  nonstochastic  radiation  pressure 
parameters,  thrusts,  polar  motion,  and  gravity  field  model  coefficients. 

The  main  requirement  on  the  reference  trajectories  is  that  they  are  within  the 
linear  remon  of  convergence  relative  to  the  true  trajectory,  since  a  linearized  (not 
extended)  Kalman  filter  has  been  adopted.  However,  if  a  reference  trajectory  is  not 
within  the  linear  region  the  Filter  can  be  used  in  a  batch  emulation  mode  to  obtain 
an  improved  set  of  initial  conditions  for  reintegration.  Reference  1  contains  a  de¬ 
tailed  description  of  the  reference  trajectory  integration  procedures  and  the  associ¬ 
ated  force  models. 

The  Corrector/Editor  system  of  programs  has  two  purposes.  The  first  purpose  is 
to  correct  tiie  data  to  instantaneous  geometric  range  (or  range  differenceljolus  clock 
offsets,  measurement  noise,  and  residual  unmodeled  and  random  effects.  Inis  is 
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accomplished  by  adjusting  the  measurements  to  account  for  time  transmission 
effects,  the  effect  of  the  displacement  of  the  antenna  system’s  electrical  phase  center 
from  the  satellite’s  center  of  gravity,  tropospheric  refraction  effects,  the  periodic 
component  of  the  relativistic  effect  on  the  satellite  clock,  and  the  solid  earth  tide 
effect  on  the  station  height.  Ran|{e  and  range  difference  data  can  be  processed  by  this 
system  of  programs  with  corrections  appropriate  to  each  type  and  source  of  data 
made.  The  second  purpose  is  to  edit  the  data.  The  Editor  does  polynomial  fits  to 
residuals,  formed  b^  differencing  a  computed  value  based  on  the  reference  trigectory 
and  nominal  clock  information  with  the  corrected  observed  value.  Consisteni^ 
checks  of  the  after-fit  residuals  are  used  to  identify  bad  points.  By-products  or  this 
editing  procedure  are  estimates  of  the  nominal  clocks  for  each  satellite  and  station, 
as  well  as  the  time  of  occurrence  and  approximate  magnitude  of  any  time  or  fre- 

3uency  Jump  events  that  may  be  observed  in  the  residuals.  Reference  2  contains  a 
etailed  description  of  the  Corrector/Editor  system  of  programs. 

Once  the  range  data  has  been  corrected  and  edited,  another  program  can  be 
used  to  search  for  simultaneous  observations  and  form  the  derived  interferometric 
data  types-differenced  range  and  doubly-differenced  range.  Differenced  range  is 
obtained  by  differencing  ranges  from  two  stations  to  the  same  satellite.  This  ^iini- 
nates  the  satellite  clock  from  the  observations.  Doubly*differenced  range  is  obtained 
by  differencing  two  differenced  ranges  involving  the  same  pair  of  stations  but  differ¬ 
ent  satellites.  This  also  eliminates  the  station  clocks  from  the  observations. 

In  addition  to  the  reference  trajectories,  edited  observations  (including  obser¬ 
vation  sigmas  and  station  coordinates),  and  nominal  clock  information  (including 
events  and  offsets  between  the  master  or  reference  clock  and  GPS  time),  the  MSF^ 
system  of  programs  requires  various  iimuts  including:  overall  program  (low  and 
identifying  information,  Quantities  defining  time  spans  for  the  fit  and  each  data 
tvpe,  lists  of  satellites  ana  stations  to  be  processed,  data  deletion  criteria,  minimum 
observation  sigmas  for  each  data  type,  the  parameter  set  and  o  priori  statistics,  and 
information  for  SATRACK  processing  if  required.  Output  products  of  the  MSF/S 
system  of  programs  include:  propagated  tregectories,  satellite  clock  offsets  from  GPS 
time  (both  time  and  frequency),  improved  initial  conditions  to  initiate  follow-on 
processing,  plots  of  corrections  and  their  corresponding  sigmas,  residual  and  signal- 
to- noise  plots,  correlation  coeflicient  matrices,  SATRACK  intersatellite  and  inter¬ 
time  covariance  matrices,  updated  station  coordinates,  and  updated  polar  motion 
information. 


ESTIMATION  CONCEPTS 


The  GPS  MSF/S  estimation  procedure  can  be  viewed  as  an  adjustment  of  a 
model  to  best  fit,  in  a  minimum  mean-squared  error  sense,  the  available  observa¬ 
tions  which  are  a  function  of  that  model.  The  possible  observations  are  range,  range 
difference,  differenced  range,  and  doubly-differenced  range-all  involving  tracking  of 
the  satellites  by  stations  on  the  surface  of  the  Earth,  The  model  consists  m  a  trsdec- 
tory  model,  clock  models  for  both  the  satellite  and  station  clocks,  and  various  ouier 
parameters  related  to  the  measurements. 

Since  the  ordinary  differential  equations  describing  a  satellite  trajectory  are 
nonlinear  and  the  Kalman  filter  equations  assume  a  linear  state  model  (and  a  linear 
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measurement  model),  a  linearisation  about  a  reference  trajectory  must  be  per¬ 
formed.  In  addition,  since  atomic  clocks  sometimes  exhibit  anomalous  step  chanma 
in  time  and/or  frequency  or  may  be  adjusted  deliberately,  a  nominal  clock  including 
approximate  values  for  all  step  changes  is  also  required,  so  the  stochastic  clock 
models  do  not  have  to  accomodate  large  jumps.  Linearisation  means  that  the  states 
of  the  Kalman  filter  are  actually  corrections  to  the  nominal  model  parameters  and 
that  the  measurements  are  processed  as  residuals.  Therefore,  partial  derivatives  are 
required  that  relate  the  states  atone  time  to  anoUier  time  (state  transition  matrices) 
and  that  relate  the  measurements  to  the  states  (measurement  matrices).  All  partial 
derivatives  are  evaluated  based  on  the  reference  trajectory  and  the  nominal  param¬ 
eter  values.  There  is  no  relineariiation  of  the  measurement  model  relative  to  the 
estimated  states.  The  A  notation  is  used  throughout  this  report  to  indicate  that 
corrections  to  nominal  model  parameters  are  actually  being  estimated. 

It  is  important  to  realise  that  a  Kalman  filter  is  a  combination  of  a  parameter 
estimation  technique  and  a  set  of  e()uations  that  define  how  the  state  ana  its  associ¬ 
ated  covariance  at  one  moment  in  time,  tj,  are  related  to  the  same  quantities  at 
another  time,  tj>|.  The  general  form  of  the  stochastic  state  equations  (in  discrete 
terms)  is  given  by: 


=  ^j^Xj  +  CVj 


(1) 


where  Axj  =  state  at 

(pj  =  =  nonsingular  transition  matrix  relating  the  state  at  tj  to 

the  state  att^+i 

Wj  =  vector  of  white  process  noise  terms  with  nonsingular  covariance 
matrix  Qy,  dim  w  ^dim  Ax 

and  C  =  matrix  of  cnes  and  zeroes  required  to  make  dim  Gw  =  dim  Ax. 

This  general  form  is  assumed  when  describing  the  processing  steps  in  a  standard 
Kalman  filter.  The  specialized  form  of  the  state  equations  adopted  in  the  MSF/S 
system,  their  underlying  models,  and  the  corresponding  Q  matrices  are  described  in 
the  STATE  EQUATIONS  section. 

The  discrete  form  of  the  linear  measurement  model  is  given  by: 

f/  =  M*v+y»  (2) 


where 

*j 

=  measurement  vector  at 

=  measurement  matrix  at  ty 

and 

=  measurement  noise  vector  at  tj. 

It  is  assumed  that  the  observations  have  been  whitened  and  decorrelated,  so  the 
measurement  noise  covariance  matrix  is  the  identity  matrix,  i.e.  =  1.  It  is  also 
assumed  that  the  processjnoise,  and  measurement  noise,  y,,  are  uncorrelated  and 
that  the  state  estimate,  Ax^,  and  its  associated  covariance  matrix,  Pq,  are  given  at  to. 
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Smoothing 


to 


Propagation 


Measurement  Update 

y/ 

{ ty}  =  mini-batch  times  ( ]  =  mini-batch  intervals  At  =  t^^;  — 


FIGIJKK  2.  SIMPlJKIKDTIMKLINK  DHPINITION 

Fi^re  2  gives  a  simplified  timeline  definition.  A  complete  timeline  definition 
is  given  in  the  next  section.  For  a  standard  linear  Kalman  niter  the  processing 
starts  by  initialisation  of  the  state  and  its  associated  covariance  matrix  at  to*  tne 
epoch  of  the  fit  span,  and  steps  forward  in  time.  Assume  for  illustration  purposes 
that  observations  only  occur  at  the  times  |  tj }.  In  the  following  equations  -  indicates  a 
predicted  quantity,  *  indicates  a  filtered  quantity,  and*  indicates  a  smoothed 
Quantity.  A  measurement  update  is  done  at  tj  by  forming  the  predicted  residuals  and 
tneir  covariance  matrix,  computing  the  gain  matrix  Kj,  and  using  these  quantities 
along  with  the  measurement  and  measurement  matrix  to  update  the  estimates  of  the 
state  and  the  covariance  matrix. 


Axj  =  Ax;  +  Kj(Xj  -  AjiiXj) 
Pj  -  a-KjAj)Pj 

where  Kj  =  Pj  AjfAjPj  Aj'h  I  )  ' 


(3) 

(4) 

(5) 


z,-A,Ax,  is  the  predicted  residual  vector  and 
AjPj  a/+  I  is  the  predicted  residual  covariance  matrix. 


The  next  step  is  to  propagate  the  state  and  the  covariance  matrix  to  time  The 
state  propagation  uses  the  state  transition  matrix  and  assumes  that  the  process 
noise  term  is  zero.  The  covariance  matrix  is  propagated  by  doing  a  deterministic 
mapping  using  the  <1>;  matrix  and  then  adding  in  Qj,  tne  state  or  process  noise  covar¬ 
iance  matrix.  This  Qj  matrix  is  the  only  difference  between  a  Kalman  filter  and  a 
sequentially  implemented  batch  weighted  least  squares  estimator. 
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The  measurement  update/propagate  pair  of  operations  is  continued  until  a  measure¬ 
ment  update  has  been  done  at  t^. 

Smoothing  is  then  accomplished  by  recursively  computing  state  ant.  covariance 
estimates  backwards,  one  step  at  a  time,  using  the  Rauch-Tung-Striebel  (RTS)  equa¬ 
tions.  Smoothed  state  estimates  at  tj  are  the  niter  estimates,  plus  a  smoother  gain 
matrix  multiplied  times  the  difference  of  the  smoothed  and  predicted  states  attj+j. 
The  gain  matrix  is  a  function  of  the  filtered  covariance  estimate  at  tj,  the  predicted 
covariance  estimate  at  tj  + ;,  and  the  state  transition  matrix  ^j.  The  smoothed 
covariance  estimate  is  the  Hltered  estimate,  plus  a  term  that  is  a  function  of  the 
difference  of  smoothed  and  predicted  covariance  estimates  at  tj+ ;  and  the  gain 
matrix  Cj. 


where 


A*; 

p; 


Axj  +  Cj  ( Ax;+ ;  -  A~Xj  +  / ) 

Pj  +  Cj‘ 


(8) 

(9) 

(10) 


This  recursive  process  is  continued  until  to  is  reached.  Reference  3  contains  an 
excellent  introduction  to  linear  Kalman  filtering  and  smoothing. 


Kalman  filters  can  easily  be  restarted  at  any  time  in  the  middle  of  the  Ht  span. 
Assuming  that  the  filter  was  previously  stopped  after  p.  '^agetion  from  t(.i  to  tt, 
restarting  the  filter  simply  consists  of  initializing  the  st,  and  covariance  matrix 
estimates  with  their  predicted  values  at  U  and  then  doing  a  measurement  update  at 
tf.  This  procedure  has  been  adopted  in  the  MSF/S  system. 


In  the  case  of  GPS  it  is  necessary  to  use  the  mini-batch  concept,  because  not  all 
observations  lie  exactly  at  the  tj  times  and  reducing  the  number  or  steps  required 
increases  the  efficiency  of  the  computations.  In  the  mini-batch  concept  all  observa¬ 
tions  in  the  interval  ^  tj  — are  processed  in  a  batch  mode.  This 


assumes  that  the  process  noise  contribution  to  state  uncertainties  can  be  ignored 
over  much  shorter  periods  than  the  time  constants  of  the  stochastic  processes.  This 
means  that  the  state  noise  covariance  matrix  Qj  is  only  added  into  the  covariance 
matrix  estimate  when  propagating  from  tj  to  ti+j.  This  essentially  averages  out  the 
random  ejects  over  the  interval  chosen,  which  primarily  affects  ^e  clock  estimates 
for  the  intervals  used  for  GPS  ( ^  1  hr).  Another  effect  is  that  solutions  are  only 
available  at  the  mini-batch  steps.  However,  that  part  of  the  state  equations 
involving  orbit-related  parameters  was  chosen  so  that  deterministic  propagation  of 
the  trajectory  corrections  between  tj  and  tj+i  is  exact. 


This  approach  also  allows  the  mini-batch  step  to  be  changed  in  the  middle  of 
the  fit  span.  This  technique  was  adopted  for  the  SATRACK  application  of  the  MSF/S 
system.  Since  this  application  requires  optimum  and  dense  estimates  only  during  a 
subspan  of  the  entire  fit  span,  a  reduced  mini-batch  step  size  span  is  defined.  The 
concepts  given  above  apply,  except  At  is  reduced  to  a  smaller  value.  The  transition 
regions  are  handled  to  ensure  no  observation  is  processed  twice.  This  technique  is 
also  useful  for  looking  in  detail  at  a  particular  time  span  to  locate  precisely  when  an 
anomaly  occurred. 
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The  square  root  information  implementation  of  the  estimation  equations  was 
selected  for  the  MSF/S  system.  This  implementation  is  mathematically  equivalent 
to  tlie  classical  Kalman  fllter/RTS  smoother  approach.  A  square  root  implementa¬ 
tion  was  selected  because  of  its  characteristic  accuracy  and  stability.  Accuracy 
involves  susceptibility  to  roundoff  errors;  stability  involves  accumulated  roundoff 
errors  not  causing  the  algorithm  to  diverge.  These  are  common  problems  when 
using  the  classical  Kalman  filter  equations  with  large  parameter  and  observation 
sets.  The  information  form  (synonymous  with  normal  equations  and  opposite  to  the 
covariance  form)  was  selected  because:  primary  interest  is  in  the  smoothed  results 
(in  this  implementation,  filter  state  and  covariance  estimates  never  need  to  be 
computed);  the  smoother  equations  require  the  inverse  of  only  an  upper  triangular 
matrix  (of  reduced  size  if  covariance  information  is  not  required),  instead  of  a  full 
matrix  (as  iti  the  RTS  formulation);  and  the  observations  are  assumed  to  already  be 
edited.  Editing  can  be  done  in  a  Kalman  filter  by  comparing  the  predicted  residual 
to  the  square  root  of  its  variance.  However,  these  quantities  are  never  computed 
explicitly  in  the  information  form.  In  fact,  no  predicted  states  and  covariances  are 
normally  computed  in  this  implementation. 

The  square  root  information  fllter/smoother  (SRIF/S)  is  based  on  the  equivalent 
concepts  of  a  data  equation  and  an  information  array  as  follows: 

Data  equation  »  Information  array 


z  -  RAx+  V  (R  z)  (11) 

where  Ax  =  states  to  be  estimated 

R  =  nonsingular  square  matrix 

V  =  zero-mean  noise  with  unity  (identity  matrix)  covariance 
z  =  right-hand  side  of  linear  equations 

Every  data  equation  corresponds  to  an  estimated  state(Ax)-covariance  (P)  pair,  i.e.. 

Ax  =  R'z  and  P  =  R'R  ’  (12) 


The  information  matrix  (normal  matrix)  is  P  —  R^R.  Therefore,  R^is  the  scmare  root 
of  the  information  matrix-the  origin  of  the  name  of  this  implementation.  Data 
equations  are  not  unique  because  if  T  is  an  orthogonal  matrix 

Tz  =  TRAx  +  Tv  (13) 

where  Tv  has  unity  covariance  and 

Ax  =  (TRjf'  Tz  -  R  't’Tz  =  R'z 
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Therefore,  the  transformed  equations  have  the  unity  covariance  noise  term  and  the 
same  solution.  These  results  can  be  extended  co  the  case  where  R  has  been  aumen- 
ted  by  additional  rows  representing  nev^  observations.  In  the  SRIF/S  method,  House¬ 
holder  orthogonal  transformations  are  used  to  partially  or  totally  trian^larize 
certain  information  arrays.  Solutions,  therefore,  are  always  computed  from  an 
upper  triangular  set  of  linear  equations.  The  Householder  transformations  are  not 
computed  explicitly  ,  only  the  results  of  applying  the  transformations  are  needed. 

The  details  of  applying  these  transformations  are  given  in  the  FILTER  ALGO¬ 
RITHM  section.  An  excellent  description  of  the  square  root  information  concepts 
and  properties  of  the  Householder  transformations,  on  which  the  MSF/S 
development  was  based,  is  in  Reference  4. 

The  parameter  set  selected  for  the  MSF/S  system  can  be  divided  into  three 
categories: 

1.  Stochastic  parameters  (labelled  p  parameters) 

a.  Orbit-related 

b.  Measurement-related 

2.  Time-varying  but  nonstochastic  parameters  (labelled  x  parameters) 

3.  Bias  parameters  (labelled  y  parameters) 

a.  Stotion-related 

b.  Orbit-related 

All  parameters  are  arc  parameters,  i.e.,  there  are  no  pass  parameters.  This  grouping 
results  in  some  simplification  of  the  SRIF/S  algorithms.  Since  all  p  parameters  must 
be  present  twice  in  the  propagation  and  array  smoothing  step  arrays,  not  allowing 
all  parameters  to  be  stochastic  results  in  smaller  arrays.  Also,  the  y  parameters  can 
be  teeated  somewhat  separately,  resulting  in  additional  array  storage  reduction.  In 
addition,  the  Q  matrix  is  only  required  for  the  p  parameters. 

Both  the  orbital  element  states  and  the  clock  states  are  treated  as  pseudoepoch 
state  parameters.  This  means  they  are  epoch  state  corrections  that  would  have 
occurred  had  the  process  noise  been  zero.  These  can  then  be  readiljr  mapped  to  cur¬ 
rent  state,  using  the  standard  state  transition  matrices.  This  definition  was  prima¬ 
rily  adopted  for  two  reasons:  (1)  to  use  the  partial  derivatives  of  position  and  velocity 
with  respect  to  orbital  elements  generated  in  the  integrator,  as  tney  would  be  used  in 
a  standard  batch  fit  when  forming  the  measurement  matrices  and  (2)  to  reduce  the 
clock  model  to  a  polynominal  in  time,  with  the  fit  epoch  as  the  reference  point  if  its 
process  noise  terms  were  zero.  This  simplifies  the  state  equations  and  keeps  the 
observation  equations  identical  to  those  used  for  batch  least  squares. 

Another  concept  employed  in  the  MSF/S  system  development  was  to  solve  for  a 
given  parameter  set  using  one  set  of  observations,  and  then,  use  these  solutions  in 
processing  another  set  of  observations.  The  two  specific  parameter  sets  chosen  were 
the  orbit-related  parameters  and  the  satellite  clock  parameters.  The  orbit  solution 
can  be  incorporated  by  using  the  propagated  trajectories  from  a  previous  fit  and  not 
solving  for  any  orbit-related  parameters.  The  satellite  clock  solutions  can  be  incor¬ 
porated  by  using  the  total  satellite  clock  offsets  from  a  previous  fit  and  not  solving 
for  satellite  clock  parameters.  These  two  methods  used  together  would  allow  sta¬ 
tions,  for  which  observations  may  or  may  not  have  been  included  in  the  previous  fit, 
to  be  positioned  with  fixed  orbits  and  satellite  clocks  as  an  evaluation  procedure. 
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TIMELINE  DEFINITIONS 


Fi^re  3  gives  a  detailed  diagram  of  the  relationships  of  all  the  time-related 
q[uantities  to  be  referred  to  in  the  rest  of  this  report.  All  times  are  either  GPS  or  UTC 
umes. 


l-f-/ — I — }■ 

to  tj  tj  +  j 


VA 


To 


■IhahI  I  -1^-//  I  ^1 . 1  I  IhHH 


Ta  te" 
FIGURES  DETAILED  TIMELINE  DEFINITION 


to 

tN 

{t/} 

At 

(] 

To 

Te 

AT 

R€d. 


=  start  time  of  nt  span 
=  end  time  of  fit  span 
=  mini-batch  times 


must  be  trajectory  timelines 


=  tj+i-tj-  mini-batch  time  step  (integer  multiple  of  trajectory 

time  step) 

=  mini-batch  measurement  limits 

=  trajectory  epoch  (may  be  different  for  each  satellite),  To^to 
=  end  time  '/f  trajectory,  Te^tjv 
=  trajectory  time  step 

=  clock  epoch  (may  be  different  for  each  satellite  and  each  station 
and  may  be  before  or  after  to) 


start  time  of  reduced  mini-batch  step  s; 
SATRACK),  =  to  +  integer  x  At 

=  end  time  of  reduced  mini-batch  step  span, 

“  ts'^***  +  integer  x  At 

=  reduced  mini-batch  time  step, 

I  integer,  At^**^  =  integer  x  AT 

=  times  for  SATRACK  covariance  information 


t^nea.  =  Start  time  of  reduced  mini-batch  step  span  (primarily  for 

^^Ued. 

{Ti}i=i;Z^ 


10 


NSWCTR  87-187 


In  addition  each  data  type  can  be  processed  over  a  subspan  of  the  fit  span  with  the 
default  being  the  entire  fit  span. 


STATE  EQUATIONS 


The  stochastic  state  equations  in  discrete  form  define  the  relationship  between 
the  states  at  tj  and  /;  i.e.,  the  solutions  that  are  determined  by  the  Filter  or 
Smoother  must  satisfy  these  equations.  The  general  form  of  these  equations  was 
^ven  in  equation  (1)  m  the  ESTIMATION  CONCEPTS  section.  A  ^ecialized 
form  of  these  equations  was  selected  for  the  MSF/S  system  to  ensure  eftlcient 
handling  of  the  multisatellite  capability  and  maximum  utilization  of  available 
partial  derivatives  and  other  pre-computed  quantities  from  the  integrator.  These 
equations  are  given  by; 

/  ^P\  /MO  0  \  /Ap  \  /  wA 

I  Ax  j  =  f  Vp^  /  0  j  J  Ax  J  -I-  I  0  J  (15) 

\Ay/  Vo  0  1/  \^y/  \o  / 

j  +  t  j 

where  p,  x,  and  y  refer  to  categories  of  parameters. 

Ap  =  stochastic  parameter  states-only  states  being  driven  by  white  noise 
Ax  =  time- varying  but  nonstochastic  parameter  states 
Ay  =  bias  parameter  states 

wj  =  white  noise  vector  with  covariance  matrix  Qj. 

The  matrix  Involving  M,  Vp ,  and  the  two  identity  matrices  is  the  state  transition 
matrix.  The  M  matrix  relates  the  p  parameter  states  at  tj  to  ;.  The  Vp  matrix 
relates  the  p  parameter  states  at  tj  to  the  x  parameter  states  at  The'^ p  param¬ 
eters  are  divided  into  two  subcategorles-orbit-related  and  measurement-related. 
The  orbit-related  p  parameters  are  radiation  pressure  (K^)  and  ^avitational  accele¬ 
rations  (G).  The  measurement-related  p  parameters  are  tropos^eric  refraction  (Cr), 
satellite  clock  (Csv)>  and  station  clock  (Cms).  The  only  x  parameters  are  the  pseudo¬ 
epoch  orbital  elements  (e).  In  addition,  the  y  parameters  are  divided  into  two 
subcategories  also-station-related  and  orbit-related.  The  only  station-related  y 
parameters  are  station  coordinates  (S).  The  orbit-related  y  parameters  are  radiation 
pressure  (RP),  thrust  (T),  pou*.  *notion  (PM),  and  gravity  coefficients  (GC).  Polar 
motion  is  under  the  orbit-related  category,  instead  of  the  station-related  category, 
because  the  partial  derivatives  of  satdlite  position  with  respect  to  polar  motion  can 
be  nonzero.  Parameters  are  present  for  either  all  satellites  or  all  stations,  depending 
on  the  specific  parameter,  except  for  thrusts  that  are  onl}'  present  for  the  appropriate 
satellites  and  polar  motion  ancl  gravity  coefficients,  which  are  common  to  all  satel¬ 
lites  and  are  present  once.  Appendix  A  describes  the  assumptions  and  definitions 
made  in  adopting  the  specialized  form  of  the  state  equations  for  the  orbit-related 
parameters. 
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The  state  transition  matrix  M  in  equation  (15)  has  the  following  block  diagonal 
structure: 


where  Nsv  =  number  of  satellites  and  N^s  =  number  of  stations. 
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The  matrix  is  also  sparse  and  has  the  following  structure,  where  each  ^  submatrix 
is  of  dimension  6X3: 
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The  Q  matrix  corresponding  to  the  p  parameters  is  also  sparse  and  has  the  following 
block  diagonal  structure: 


The  Filter  algorithm  does  not  directly  use  the  Q  matrix;  but  rather,  uses  the  square 

root  of  the  inverse  of  the  Q  matrix  i.e.,  Q  =  R'l  rZ.  Rw  bas  the  same  block 

diagonal  structure  as  Q  except  that  each  block  is  upper  triangular. 

The  parameters  are  detailed  below.  It  is  assumed  that  the  values  solved  for  are 
corrections  to  the  nominal  values  for  all  parameters.  Included  in  the  description  of 
each  parameter  are  the  models  assumed  and  a  definition  of  the  submatrices  of  M,  Vp, 
and  Q,  if  applicable.  Also,  units  are  nven  for  each  parameter.  The  x  parameters  are 
described  nrst,  so  Uie  submatrices  ofVp  can  be  defined  for  the  p  parameters  as  they 
are  described.  Several  expressions  involve  At  and  must  be  re-evaluated  for  the 
reduced  mini-batch  span  processing. 


14 


NSWCTR  87-187 


TIME-VARYING  NONSTOCHASTIC  STATES 

The  only  x  parameters  are  the  pseudoepoch  orbital  elements  at  the  trejectory 
epoch  represented  by  a.  Each  trajectory  can  nave  a  different  epoch;  as  long  as  it  is 
before  the  fit  epoch.  These  orhital  elements  are  defined  as  follows: 

a  =  semi-ms^or  axis  (km) 

e  sin  (0  =  eccentricity  X  sine  (argument  of  perigee)  (unitless) 
e  cos  u  =  eccentricity  X  cosine  (argument  of  perigee)  (uuitless) 

I  =  inclination  (radians) 

M  +  (o  =  mean  anomaly  +  argument  of  perigee  (radians) 

Q  =  right  ascension  of  the  ascending  node  (radians). 

The  inertial  position  ro  and  velocity  fo  can  alternatively  be  the  x  parameters,  if 
partials  with  respect  to  these  parameters  are  on  the  tredectory  for  any  given  satellite. 
The  MSF/S  system  of  programs  is  designed  so  that  it  does  not  matter  which  of  these 
two  parameter  sets  is  chosen  for  a  given  satellite.  Throughout  the  rest  of  this  report 
o  will  be  used  to  signify  either  set  of  orbit  parameters.  No  direct  process  noise  is 
included  on  the  a  parameters,  nevertheless,  they  are  smoothable  as  a  result  of  being 
dynamically  related  (through  the  Vp  matrix)  to  the  Kr  and  G  parameters,  which  are 
stochastic.  Because  of  the  pseudoepoch  state  formulation  of  the  state  equations  (see 
Appendix  A),  the  state  transition  matrix  for  these  states  is  the  identity  matrix. 


STOCHASTIC  STATES 

As  mentioned  above,  the  p  parameters  can  be  divided  into  two  subcategories- 
orbit-related  and  measurement-related.  Since  all  x  parameters  are  orbit  param¬ 
eters,  ordering  the  orbit-related  p  parameters  first  results  in  the  Vp  matrix  being  of 
the  form  given  in  equation  (17). 


Orbit-Related  Stochastic  States 


Radiation  Pressure 


The  3  radiation  pressure  parameters,  Kr  ,  are 

/  radiation  pressure  scale  (unitless,  .01  constant  acceleration  \ 


/ 

Kg, 

= 

\  K«J 

of  10'*^  km/sec^  in  satellite-sun  direction) 
y-axis  acceleration  (10  km/sec^ ) 
angle  between  satellite  x  and  y  axes  (radians,  usually  constrained 


i  angle  between  satellite  x  and  y  axes  iradi 
\  to  90  deg  since  it  is  not  easily  observed). 


/ 


15 


NSWCTR  87-187 


Nominal  values  for  these  parameters  are  those  used  in  the  intenator.  The  time 
history  of  the  corrections  to  each  of  these  parameters  is  modeled  as  a  current  state 
first-order  Gauss-Markov  process  (see  Appendix  B).  This  results  in  the  3X3  portion 
of  the  state  transition  matrix  M  being  given  by 


(19) 


where  e~  reduces  to  1  for  =  *•  i-e.,  for  a  random  walk  process  or  bias.  The 
process  noise  covariance  matrix  for  these  states  is  a  diagonal  matrix  with  each 
element  given  by 


i  =  1,2,3 


(20) 


A 

where  =  spectral  density  of  the  white  noise  term  =  ^ and  is  the 

steady-state  sigma  of  the  process.  These  si^as  have  the  same  units  as  the  param¬ 
eters.  Each  dia^nal  of  the  matrix  for  these  states  is  the  square  root  of  the 
inverse  of  each  diagonal  of  the  Qic„  matrix.  For  a  random  walk  process  the  spectral 
density  must  be  directly  input.  This  element  of  the  process  noise  covariance  matrix 
expression  ’‘educes  to  At  for  a  random  walk  process  and  0  for  a  bias,  in  which  case 

the  diagonal  of  Ru,  is  set  io  a  computational  infinity.  The  6X3  portion  of  the  Vp 
matrix  corresponding  to  Kk  is  given  by 


/  dritj^i) 
f  aK«(tj) 
I  ar(tj^;) 


(21) 


where  4>e  (tj+ ;)  —  partials  of  position  and  velocity  at  time  t^ + ;  with  respect  to  epoch 
orbital  elements  obtained  by  interpolation  from  the  trajectory. 


Approximations  using  second  order  Taylor  series  expansions  are  used  to  obtain  the 
partials  of  position  and  velocity  at  time  with  respect  to  radiation  pressure 
parameters  at  tj.  (This  approximation  is  discussed  in  Appendix  C.)  These  partials 
are  given  by 


2  dKR(tj) 


(22) 


akfiTt^T 


ar^t.)  At*  apft,) 

sKrUj)  2  dK^dj) 


(23) 
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where 


0 


0 


(24) 


=  0  for  a  random  walk  or  bias 


j  ag-Kn^IO  shape  cos  10 shape  cos  -Kp^I0'‘*  shape  sin  Kp^ 

Ki^, 


dritj) 

dKpitj) 


0 

— 


10 shape  sin  Kp^  KpJ0'‘^  shape  cos  Kp^  j  (26) 
0  0 


a  =  inertial  acceleration  at  tj  due  to  radiation  pressure  in  the  body-axes 
directions  obtained  from  the  trigectory 

R,  -  matrix  required  to  transform  between  the  body-axis  and  inertial 
Cartesian  reference  systems  obtained  from  the  trigectory  at  tj 

Kp  -  nominal  radiation  pressure  parameter  values  from  the  trigectory 

shape  =  fraction  of  the  sun's  disk  unobstructed  by  any  eclipsing  body  (Earth, 
Moon,  or  both)  obtained  from  the  trcgectory  at  tj 

The  effects  of  changes  in  e  on  Kp  are  ignored  because  they  are  negligible.  Kp  parame¬ 
ters  cannot  be  present  without  e  parameters.  The  partials  of  radiation  pressure 
acceleration  wiUi  respect  to  Kp  will  be  0  if  lies  in  the  umbra  region  of  the  eclipse, 
even  though  some  observations  in  the  corresponding  mini-batch  interval  may  lie 
outside  of  the  eclipse  period.  This  form  of  mismodelmg  is  not  a  problem  when  using 
these  parameters  in  tneir  intended  stochastic  manner. 


Gravitational  Accelerations 


The  3  gravitational  acceleration  parameters,  G,  are  the  radial,  along- track,  and 
cross-track  O^C)  accelerations  in  km/sec^.  They  are  called  this  because  we  primary 
error  they  are  meant  to  absorb  is  cavity  Held  model  error.  These  parameters  are 
sometimes  referred  to  as  unmodeled  aixelerations.  The  nominal  value  of  each  of 
these  parameters  is  assumed  to  be  zero.  The  time  history  of  the  corrections  to  each  of 
these  parameters  is  modeled  as  a  current  state  first-order  Gauss-Markov  process  (see 
Appendix  B).  This  results  in  the  3X3  portion  of  the  state  transition  matrix  M  being 
given  by 
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Mg  = 


(26) 


At 

where  e  "  reduces  to  1  for  xc  =  *,  i.e.,  for  a  random  walk  process  or  bias.  Th  j 
process  noise  covariance  matrix  for  these  states  is  a  diagonal  matrix  with  each 
elem«int  given  by 


/  2At\ 

Qc,,  =  Qg,  "■  *  y  *  “ 

2 

where  qc,  =  spectral  density  of  the  white  noise  term  =  -  -  Og,  and  o^,  is  the  steady- 

state  sigma  of  the  process.  These  sigmas  have  the  same  units  as  the  parameters. 
Each  diagonal  of  the  Rm  matrix  for  these  states  is  the  square  root  of  the  inverse  of 
each  diamnal  of  the  Qg  matrix.  For  a  random  walk  process  the  spectral  density  must 
^  directly  input.  The  process  noise  covariance  matrix  expression  reduces  to  qc  At 
for  a  random  walk  process  and  0  for  a  bias,  in  which  case  uie  diagonal  of  R„  is  set  to  a 
computational  infinity.  The  6X3  portion  of  the  Vp  matrix  corresponding  to  G  is  given 
by 


dCitj)  \ 

j  (28) 

dCitj)  J 

where  <(>,  (tj+j)  =  partials  of  position  and  velocity  at  time  tj^-j  with  respect  to  epoch 
orbital  elements  obtained  by  interpolating  off  the  trajectory. 

Approximations  using  second  order  Taylor  series  expansions  are  used  to  obtain  the 
partials  of  position  and  velocity  at  time  tj+i  with  respect  to  gravitational  accelera¬ 
tion  parameters  at  tj.  (This  approximation  is  discus^  in  Appendix  C.)  These 
partials  are  given  by 


_  At*  dir(/,) 

(29) 

dCitj) 

2  dC(tj) 

a?  ft,)  At*  ^(tj) 

“At  +  Bn 

(30) 

dG(tj) 

aery  2  dc(tj)  ® 
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where 


(31) 


=  0  for  a  random  walk  or  bias 

r  <32) 

acfi/)  V  iFx^  ir^/ 


=  matrix  required  to  transform  between  the  RAC  and 
inertial  Cartesian  reference  frames  at  tj,  where  * 
denotes  a  unit  vector 

The  effects  of  chanms  in  a  on  G  are  ignored  because  they  are  negligible.  C 
parameters  cannot  m  present  without  e  parameters. 


Measurement-Related  Stochastic  States 
Tropospheric  Refraction 

The  tropospheric  refraction  parameter,  Cm,  is  the  xenith  tropospheric  refraction 
parameter  in  km.  The  correction  to  this  parameter  is  related  to  other  elevations  by 
a  factor  of  1/dn  E.  The  time  history  of  the  correction  to  this  parameter  is  modeled  as 
a  current  state  first-order  Gauss-Markov  process  (see  Ai>pendix  B).  This  results  in 
the  1 X 1  portion  of  the  state  transition  matrix  M  being  pven  by 


At 


Ale,  =  •  ■  xc^ 


(33) 


where  •  "  xc,  reduces  to  1  for  xc=  i.e.,  for  a  random  walk  process  or  bias.  The 
process  noise  covariance  matrixTor  this  state  consists  of  a  single  element  given  by 


Qc,  -  Qc,  ^ 


(34) 


2 

where  qc,  =  spectral  density  of  the  white  noise  term  3=  —  and  oc^  is  the 

steady-state  siraa  of  the  process.  This  sigma  has  the  same  unit  as  the  parameter'. 
The  Ru,  matrix  Tor  this  state  is  the  square  root  of  tlte  inverse  of  Qc,-  For  a  random 
walk  process  the  spectral  density  must  be  directly  input.  The  process  noise  covari¬ 
ance  matrix  expression  reduces  to  qcj,At  for  a  random  walk  process  and  0  for  a  bias, 
in  which  case  R^  is  set  to  a  computational  infinity. 
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Satellite  and  StoUon  Ciockt 

The  clock  models  for  the  satellite  and  station  clocks  are  given  below. 

Appendix  D  discusses  the  models  for  the  satellite  and  station  clocks  in  more  detail 
than  given  here  rnd  describes  the  relationriite  between  the  clock  model  spectral 
density  noise  terms  and  the  Allan  variance.  The  firequency  offset  term  in  these 
models  is  not  the  instantaneous  frequency  oflbet»  since  it  does  not  contain  the  white 
frequency  noise  component  This  component  is  only  observable  bv  its  integrated 
effect  on  the  time  offset  These  clock  models  reduce  to  polynomials  referenced  to  the 
fit  epoch  when  all  the  process  noise  terms  are  lero. 

The  satellite  clock  parameters,  Csv*  are 


/Aio\  /frequency  drift  (ppm/sec)  \ 

Csv  =  dio  =1  frequency  offset  (ppm)  I 


^Ato  j  \  time  oflset  (psec) 


The  nominal  values  for  these  parameters  are  described  in  the  next  section  and  are 
based  on  initial  pol^omials  and  step  changes.  The  clock  model  is  implemented  in  a 
pseudoepoch  state  term  which  results  in  the  portion  of  the  state  transition  matrix 
corresponding  to  these  states  being  an  identity  matrix,  i.e., 


"Csv 


(35) 


The  state  noise  covariance  matrix  for  this  set  of  states  is  given  by 

At*  \ 

Qi-=-  \ 


1  Qt  lit 

At* 

Qiy 

"Csv 

At* 

q-T 

dt* 

Qjy  +Qadt 

At" 

At*  At*  At* 

Qiy  +  q.  y 

/  ' 

where 

♦Csv=  1  + 

At^  At* 

q,—  +  q* 

%  2 


At* 

T 


•T  -I  'T 

^^sv  ”■  ^Cgv^Csv^^Csv  (38) 


0  \ 
0 


1 


) 


(37) 
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{qj}  =  white  noise  spectral  density  values  in  units  of  (ppm/sec)'*  /sec, 
(ppm)^/sec,  and  (psec)^  /sec;  and  tj+i  -to  is  in  seconds. 

The  Filter  actually  requires  To  get  this  matrix  is  factored  into 

^Csv  ^Csv  using  lower  triangular'Cholesky  decomposition  (see  Appendix  E).  Then, 
1>Csv  is  no  longer  upper  triangular  and  must  be  upper  trian^larized  before 
being  used  in  the  propagation  step  for  computational  sWplifications.  If  some  of  the 
Qi  values  are  zero,  special  processing  must  be  done  to  ensure  that  is  nonsingular. 
Ifqj  *  0  and  qs,  qj  or  both  =  0,  no  change  is  necessary.  If  q;  ~  0,  Qi, i  is  set  to  =»0.  Ifq/ 
=  Qa  =  p.  Q2JI  is  set  to  =0.  If  q  /  =  qa  =  Qj  =  0,  Qaj  is  set  to  **0.  If  only  processing 
range  difference  data,  each  is  set  to  ^0.  However,  the  a  priori  sigmas  on  the 
parameters  should  not  be  set  to  ^  0.  Because  the  time  offset  state  is  basically  the 
integral  of  the  frequency  offset  state,  constraining  the  time  offset  also  constrains  the 
frequency  offset,  which  is  not  desired. 

For  certain  clock  events  the  matrix  is  changed  for  the  mini-batch  step 

propagation  that  contains  the  event  and  then  changed  back  to  its  orimnal  values. 
This  is  used  to  account  for  the  uncertainty  in  the  clock  event  input  oi^sets. 

For  a  C-field  adjust,  R'c^y  replaces  Rcsv. where 


=z  Or^ 

^SV 
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and 


QCsv  =  QCsy  + 


where  each  input  sigma  is  converted  to  the  proper  units  before  use. 
The  station  clock  parameters,  Cms,  are 


Cms 


f  frequency  offset  (ppm)  \ 
time  offset  (psec)  j 


(42) 


This  model  is  identical  to  the  satellite  clock  model,  except  the  frecjuency  drift  state  is 
absent.  The  nominal  values  are  also  based  on  an  initial  polynomial  and  step 
changes.  The  model  is  implemented  in  a  pseudoepoch  state  form,  which  results  in 
the  portion  of  state  transition  matrix  corresponding  to  these  states  being  an  identity 
matrix,  i.e.. 


^c«c  - 


1  0 
0  1 


The  state  noise  covariance  matrix  for  this  set  of  states  is  given  by 


MS 


/  \ 

/  q/At  qi  —  \ 


At^  At^ 

Qi —  Qi —  +  q^At 
2  3 


(43) 


(44) 


where 


1 


-  to 


(45) 


{qd  =  white  noise  spectral  densities  in  units  of  (ppm)Vsec  and  (psec)Vsec;  and 
tj +/  -  to  is  in  seconds 

The  Filter  actually  requires  =  Qc^s  Ket  this  matrix  Qc^s  is  first 

MA 

7' 

factored  into  Rc„,  using  lower  triangular  Cholesky  decomposition  (see 
Appendix  E).  Then,  no  longer  upper  triangular  ana  must  be  upper 

tnangpilarized  before  being  used  in  the  propagation  step  for  computational  simpli¬ 
fications.  If  some  of  the  qi  values  are  zero,  special  processing  must  be  done  to  ensure 
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that  Qc^fs  is  nonsin^lar.  If  q/  =  0,  Qtj  is  set  to  =»0.  If  in  addition  q^  =  0,  is  set  to 
»»  0.  If  only  processing  range  difference  data,  each  q^  is  set  to  *»  0.  However,  the 
a  priori  sigmas  on  the  parameters  should  not  be  set  to  ^0.  This  is  because  the  time 
onset  state  is  basically  the  integral  of  the  frequency  offset  state.  Constraining  the 
time  offset  also  constrains  the  frequency  offset,  which  is  not  desired.  For  the  selected 
master  station  the  qi  values  are  set  to  ^0. 


For  certain  clock  events,  the  is  changed  for  the  mini-batch  step 

propagation  that  contains  the  event  and  then  changed  back  to  its  orimnal  values. 
This  is  used  to  account  for  the  uncertainty  in  the  clock  event  input  onsets.  For  a 
clock  reinitialization, replaces  Rc^g,  where 


I' 


0  \ 
J _ 


(46) 


and  each  sigma  is  the  a  priori  sigma  used  in  the  Filter  initialization.  This  is  not  done 
for  a  master  station  clock  reinitialization,  which  is  intended  to  accomodate  GPS  time 
steering.  For  a  frequency  change  event,  Rci^^  replaces  Rct,s>  where 


(47) 


and 


(48) 


where  each  input  sigma  is  converted  to  the  proper  units  before  use.  For  a  master  sta¬ 
tion  switch  event,  the  Rc„«  matrix  for  the  station  that  is  no  longer  the  master  station 
is  replaced  by 


(49) 


for  the  first  propagation  step  for  which  tj+j  >  t^ss-  Then,  the  matrix  is  reset  to  the 
Rc^  values  ori^nally  computed  from  input  and  saved  (before  it  was  replaced  by  a 
ma&ix  based  on  q/s  set  to  ”«*0)  for  propagation  from  tj+i  to  tj+2  and  all  subsequent 
steps.  The  master  station  Qc^s  (essentially  a  null  matrix)  is  then  used  from 

tj  to  tj+t  and  all  subsequent  propagation  steps  for  the  new  master  station  parame- 
teis.  No  change  in  the  process  noise  is  required  for  a  station  time  change  event,  since 
it  is  assumed  to  be  exact. 
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BIAS  STATES 

As  previously  mentioned,  the  y  parameters  are  divided  into  two  subcategories- 
station-related  and  orbit-related.  Alforbit- related  y  parameters  require  partial 
derivatives  from  the  trajectories  except  polar  motion,  in  which  case  partials  are 
present  only  if  the  geopotential  expansion  axis  is  not  the  instantaneous  spin  axis  or 
the  Celestial  Ephemeris  Pole. 

Station- Related  Bias  States 


Station  Coordinates 

The  corrections  to  station  coordinates,  AS,  are  defined  in  a  local-vertical 
reference  frame  as  follows: 

/  AE  \  /  east  component  \ 

AS  =  (  An  j  =  (  north  component  )  (km) 

\  A V  /  \  vertical  component  / 

The  nominal  station  coordinates  are  input  in  terms  of  longitude,  latitude,  and  height 
referenced  to  a  specified  ellipsoid.  The  state  transition  matrix  for  these  states  is  the 
identity  matrix. 

Orbit- Related  Bias  States 


This  category  of  parameter  is  included  so  constant  force  model  parameters 
affecting  the  orbit  can  be  solved  for.  All  of  these  parameter  sets  have  identity  state 
transition  matrices;  and  because  of  the  pseudoepoch  orbital  element  state  deHnition, 
they  result  in  no  terms  that  relate  changes  in  e  to  changes  in  these  parameters.  Par¬ 
tial  derivatives  are  only  required  in  forming  the  measurement  condition  equations 
and  are  obtained  by  interpolating  off  the  trsgectory.  Four  parameter  sets  fall  in  this 
category. 

Radiation  Pressure 


The  3  radiation  pressure  parameters,  RP,  are 


!  RP,  ^  j  radiation  pressure  scale  (unitless) 

y-axis  acceleration  (10  km/sec'O 


RPs 


angle  between  the  x  and  y  axes  (radians  are  usuallv 
constrained  to  90  deg,  since  it  is  not  easily  observed) 


/ 


These  parameters  are  present  for  every  satellite,  if  selected,  and  are  identical  to  the 
stochastic  radiation  pressure  parameters  except  the  corrections  are  modeled  as  epoch 
state  constants.  The  stochastic  parameter  states  can  be  confi^red  to  emulate  cur¬ 
rent  state  constants.  However,  this  method  is  not  reconunended  because  of  the 
partial  derivative  approximations  made  and  the  fact  that  the  smoothing  j^ocedure 
would  have  to  be  completed  to  obtain  the  orbit  solution  at  each  timeline.  The  RP 
parameters  should  be  used  in  two  cases;  to  emulate  a  batch  orbit  fit  using  the  Filter 
or  to  estimate  constants  along  with  their  stochastic  counterparts  (since  the  latter 
states  are  assumed  to  be  of  zero  mean). 
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Thrust 

The  thrust  parameters,  T,  are 


(first  component  \ 
second  component  I 
third  component  / 


(km/sec“) 


Thrusts  can  be  modeled  in  the  inte^ator  in  one  of  three  reference  frames-the  body- 
axis  frame,  the  RAC  frame,  or  the  RVC  frame  (see  Reference  1).  Thrust  parameters 
are  only  present  for  specific  satellites,  as  required,  and  their  nominal  values  are 
obtained  from  the  trajectories. 

Polar  Motion 

The  polar  motion  parameters,  PM,  are 


/ 


\ 


P 

Q 


/  pole  component  along  Greenwich  meridian  (radians) 


/ 


pole  component  along  meridian  90  deg  west  of 
Greenwich  (radians) 

\  rate  ofehange  of  UTl-UTC  (sec/sec) 


The  corrections  (o  these  par^eters  are  modeled  as  constants  over  the  entire  fit  span 
with  the  At  term  reference  time  being  the  fit  epoch.  These  states  are  common  to  all 
satellites  and  stations. 


Gravity  Coefficients 

The  gravity  coefficient  parameters,  GC,  are  selected  gravity  field  model 
coefficients  (unitless).  These  parameters  are  common  to  all  satellites  and  the 
nominal  values  are  those  used  by  the  integrator. 

When  the  orbit-related  parameters  are  viewed  together,  several  conunents 
apply.  With  only  one  set  of  orbit-related  stochastic  parameter  states  and  pseudo¬ 
epoch  orbital  element  states  being  solved  for,  the  orbit  model  is  essentially  equiva¬ 
lent  to  solving  for  position,  velocity,  and  acceleration  corrections  with  the  accelera¬ 
tion  corrections  constrained  to  be  a  zero-mean  Gauss- Markov  or  random  walk  proc¬ 
ess.  The  Ka  and  C  parameters  should  not  be  used  as  stochastic  parameters  simulta¬ 
neously.  If  K/i  is  chosen,  it  can  be  viewed  as  solving  for  acceleration  corrections  along 
the  satellite-sun  line  and  along  the  y  axis.  This  is  because  the  direct  radiation  pres¬ 
sure  force  is  almost  constant  throughout  the  orbit,  so  the  scale  parameter  is  simply 
scaling  this  near-constant  acceleration.  If  G  par  ameters  are  chosen,  the  acceleration 
solved  for  is  resolved  into  the  RAC  coordinate  frame.  The  constant  RP  parameters 
could  be  used  with  the  stochastic  Km  parameters  to  account  for  non-zero  mean  accel¬ 
erations  in  the  satellite-sun  and  y-axis  directions.  A  continuous  thrust  would  serve 
the  same  purpose  for  the  G  parameters. 
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OBSERVATION  EQUATIONS  AND  PARTIAL  DERIVATIVES 


RANGE 


As  mentioned  in  the  GPS  DATA  FLOW  section  it  is  assumed  that  all  data  to 
be  used  in  the  MSF/S  system  have  been  corrected  for  time  transmission,  relativity, 
satellite  antenna  offset,  tropospheric  refraction,  and  solid  earth  tide  station  height 
effects.  Therefore  the  observation  model  is  considerably  simplified.  For  range  (R) 
data  the  nonlinear  observation  equation  is  given  by 


Ri,*  (to6»^  ~  I  JO*  f  (to6i()J 

+  + A£»(u.)]+  sinEjltbs) 

where  i  =  satellite  subscript 
k  =  station  subscript 
tob»  ~  observation  time 

riitobs)  —  inertial  coordinates  of  satellite  t  at  tabu  interpolated  off  the 
appropriate  trajectory 


(tofts)  ~  [ABCD  (toft»)f  T/kj 

=  inertial  coordinates  of  station  k  at  tobu 


(50) 


(51) 


r* 


KF 


/  COS  cos  \k\  /  ®  \ 

=  (A'  +  h*)  f  COS sin j  —  Ae^i  0  ] 

\  sin^k  / 


=  Cartesian  Earth-fixed  station  coordinates 
(k/^,  hk)  =  geodetic  coordinates  of  station  k 

j  =  oblateness  of  Earth 

e  =  K^-nn* 


A  = 


°Eartk 


(I  —  e^sin^$*)^ 


a 


Earth 


=  semimajor  axis  of  Earth's  reference  ellipsoid 


(52) 


(53) 

(54) 


ABCD  (tab J  =  inertial-to-Earth-fixed  rotation  matrix  at  toftg.  (Reference  1 
contains  the  details  of  the  inertial-to-Barth-nxed 
transformation  computations.) 
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c  =  speed  of  light  (km/sec) 

A(f(tu6«)  —  nominal  clock  time  offset  of  the  ith  satellite  clock  from 

GPS  time  at  to6«  (to  be  described  below) 

(tabs)  =  time  offset  correction  to  the  nominal  clock  for  the  ith  satellite  at  to6« 
(nominally  zero) 

^ikitobt)  =  nominal  clock  time  offset  of  the  ftth  Station  clock  from 
GPS  time  at  t^bt  (to  be  described  below) 

(tu6<>)  =  time  offset  correction  to  the  nominal  clock  for  the  ^th  station  at  tob$ 

(nominally  zero) 

dobs)  -  zenith  tropospheric  refraction  correction  at  tabs  (nominally  zero) 
Ei.*(t(ifcs)  =  90"' —  Arccos(p  uv)  (65) 

=  instantaneous  elevation  angle  from  station  k  to  satellite  i 
P  =  -  rkdohs)  (56) 

Uv  =  [ABCD(uJV  (  yK  )  (57) 

\zef  / 

/xkF\ 

The  vector  f  ykp  1  is  computed  by  evaluating  equation  (52)  with  hk  =  0,  and 

\ZEF/k 

'  indicates  a  unit  vector. 

Computation  of  A{"(t)  and  Alf(t) 


The  nominal  satellite  time  and  frequency  offsets  at  an  arbitrary  time  t  are 
given  by 


Ati(t)  —  AtOj  -f-  Aio,  (t  —tc)  +  Alt). 


(t-tc/ 


M-it)  =  Aio.-I-  AVo,(f  -tc.)  (59) 

where  Ato ,  Aioj,  and  AVo,  are  input  quantities  converted  to  internal  units  (n8ec-»psp>.:^ 
parts  in  10'^-»ppm,  and  parts  in  lO'Vday-^pom/sec)  and  tc,  is  the  ith  satellite's  clock 
epoch.  To  accomodate  jumps  in  the  nominal  clocks,  four  satellite  clock  events  have 
been  defined.  These  events  are  described  below  in  terms  of  their  effects  on  the  nomi¬ 
nal  clocks.  The  corresponding  clock  process  noise  adjustments  were  discussed  in  the 
previous  section.  All  events  are  input  as  time  of  the  event  (in  day  number  and  sec¬ 
onds  of  the  day)  and  associated  clock  offsets  (in  nsec,  parts  in  10'^  and  parts  in 
lO'Vday).  These  offsets  are  converted  to  the  above  internal  units  before  use. 
Processing  of  each  event  is  described  below. 
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1.  O-fleld  ac^ust 

A  C-fleld  adjust  is  a  generic  name  for  a  conunanded  change  in  the  satellite's 
clock  frequency  (usually  only  applies  to  Cesium  clocks).  The  amount  of  the  change  is 
not  known  exactly,  so  a  frequency  uncertain^  is  added  to  the  process  noise  matrix 
for  this  event  (see  the  STATE  EQUATIONS  section^  This  event  affects  the 
nominal  clock  as  follows: 

For  the  first  time  t  >tc./uw, 

Aio,  is  replaced  by  -tc)  + 

is 

and  AJ(,.  is  replaced  by  Alf,.  +  +  AVo.(tc./«w.  -tc.) 

All  subsequent  nominal  clock  computations  for  this  satellite  are  referenced  to  a 
redefined  epoch,  tc  /ieu,,  and  use  equations  (58)  and  (59). 

2.  Z-count  adjust 

A  Z-count  adjust  is  a  commanded  change  in  the  satellite  clock's  time  offset.  The 
amount  of  the  change  is  known  exactly.  This  event  affects  the  nominal  clock  as 
follows: 

For  the  first  time  t  &  tz-coum^ 

At(,^  is  replaced  by  Ai(,.  +  ^iz-count, 

All  further  nominal  clock  computations  for  this  satellite  are  based  on  equations  (58) 
and  (59)  with  the  clock  epoch  unchanged. 

3.  Clock  reinitialization 

A  clock  reinitialization  is  either  a  switch  in  the  operational  clock  on  the 
Si .  Hite  or  an  anomalous  phase  jump  in  the  current  clock.  The  uncertainties  in  the 
r  1  c  estimates  are  set  back  to  approximately  the  initialization  values  (see  STATE 
E  J  Al'IONS  section).  This  event  affects  the  nominal  clock  as  follows: 

For  the  first  time  t  s  trtmu, 

Aio-,  A'  and  AVo,  are  replaced  by  a  new  set  of  values  from  input  and  the  clock  epoch 
is  redc./ned  to  trtinu,  for  all  further  nominal  clock  computations. 

Frequency  change 

A  frequency  change  is  an  unexplained  jump  in  frequency.  For  this  event,  it  is 
possible  to  increase  ^e  uncertainty  in  each  clock  state  (see  STATE  EQUATIONS 
section).  This  event  affects  the  nominal  clock  as  follows: 

For  the  first  time  t  > 

Aio-  is  replaced  by  Aio,  +  Aio^(t/'c_  -^1  +  A’io- 

£ 

At^i  is  replaced  by  AiYc, 
and  AVo,  is  replaced  by  AV/i, 
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All  subsequent  nominal  clock  computations  for  this  satellite  are  referenced  to  a 
redefined  epoch,  t,  ,  and  use  equations  (58)  and  (59). 

'Cl 

If  the  total  satellite  clock  offsets  from  a  previous  Filter  or  Smoother  run  are 
available, 

Atr(to6s)  +  Ati  (tohj  is  replaced  by  +  ^iToiai,itj){toh,-tj) 

where  ^j-  -^  <  ^obs  ^  tj  +  -^ 

The  total  satellite  clock  offsets  are  defined  in  the  SOLUTION  AND 
DIAGNOSTICS  secUon. 

Computation  of Att(t)  and  Au(t) 

The  nominal  station  time  and  frequency  offsets  at  an  arbitrary  time  t  are  given 
A(t(t)  =  Aioj  +  £iloi,(t—tcJ  +  Aioj  — -  (60) 

A 

AjJ(t)  =  AJo*  +  AV<,*rt-te,>  (61) 

where  Ato^,  Ai^j,  and  AVo.  are  input  quantities  converted  to  interna)  units 
(nsec-»usec,  parts  in  lO'^ppm,  and  parts  in  10'Vday-*ppm/sec)  and  tc^  is  the  kth 
station's  cIock  epoch.  To  accomodate  jumps  in  the  nominal  clocks,  three  station  clock 
events  have  been  defined.  These  events  are  described  below  in  terms  of  their  effects 
on  the  nominal  clocks.  The  corresponding  clock  process  noise  adjustments  were 
discussed  in  the  previous  section.  All  events  are  input  as  time  of  the  event  (in  day 
number  and  seconds  of  the  day)  and  associated  clock  offsets  (in  nsec,  parts  in  10‘S 
and  parts  in  lO'Vday).  These  offsets  are  converted  to  the  internal  units  before  use. 
Processing  of  each  event  is  described  below. 

1.  Clock  reinitialization 

A  clock  reinitialization  is  either  a  switch  in  the  operational  clock  at  the  station 
or  an  anomalous  phase  jump  in  the  current  clock.  The  uncertainties  in  the  clock 
estimates  are  setWck  to  approximately  their  initialization  values  (see  STATE 
EQUATIONS  section).  This  event  affects  the  nominal  clock  as  follows: 

For  the  first  time  t  ^ 

Aio*,  At^,  and  AVo^  are  replaced  by  a  new  set  of  values  from  input  and  the  clock 
epoch  is  redefined  to  for  all  further  nominal  clock  computations.  This  event  is 
also  used  to  accommodate  the  GPS  time  steering  procedure  when  the  offsets  between 
GPS  time  and  the  master  clock  time  change,  but  the  designated  master  clock  does 
not. 


2.  Frequency  change 

A  frequency  change  is  an  unexplained  jump  in  frequen^.  For  this  event,  it  is 
possible  to  increase  the  uncertainty  in  each  clock  state  (see  STATE  EQUATIONS 
section).  This  event  affects  the  nominal  clock  as  follows: 
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For  the  first  time  t  >  V, 

.  .. 

Aioj  is  replaced  by  Ajo*  +  Ato*  +  AVo^  — - 

Aro|^  is  replaced  by  Ai^^ 
and  Aro||  is  replaced  by  AV^^ 

All  subsequent  nominal  clock  computations  for  this  station  are  referenced  to  a 
redefined  epoch, ,  and  use  equations  (60)  and  (61). 

3.  Time  change 

A  time  change  is  a  commanded  chaiijK  in  the  station  clock's  time  offset.  The 
amount  of  Uie  change  is  known  exactly.  Tnis  event  aflects  the  nominal  clock  as 
follows: 


For  the  first  time  tztu 


Afoj^  is  replaced  by  A^o^^  + 


All  subsequent  nominal  clock  computations  for  this  satellite  are  based  on  equations 
(60)  and  (ol)  with  the  clock  epoch  unchanged. 


4.  Master  station  switch 

The  master  station  switch  event  involves  two  stations  simultaneously:  the  original 
and  the  new  master  stations.  For  the  station  that  is  no  longer  the  master,  a  clock 
reinitialization  is  dene  at  t^ss,  i.e.. 


Aro^i  Arojk,  and  AVo^^  are  replaced  by  new  input  quantities, 

the  clock  epoch  is  redefined  to  be  tmsst  and  all  subsequent  nominal  clock  computa¬ 
tions  are  referenced  to  this  redefined  epoch.  For  the  new  master  station 

Aioj^.,  Atoj^.,  and  AVo^^.  are  replaced  by  new  input  quantities, 

the  clock  epoch  is  redefined  to  be  tMss>  nnd  all  subsequent  nominal  clock  compu¬ 
tations  for  this  station  are  referenced  to  this  redefined  epoch.  The  process  noise 
terms  for  these  stations’  clocks  are  adjusted  as  specified  in  the  previous  section.  In 
addition,  the  information  array  elements  corresponding^to  the  clock  states  for  the 
new  master  station  have  to  be  modified.  (See  the  FILTER  ALGORITHM  section.) 


The  linear  measurement  model,  equation  (2),  requires  computation  of  the 
measurement  residuals,  x,  and  partial  derivatives  of  the  observations  with  inspect 
to  the  parameters,  A.  Each  residual  is  computed  by  evaluating  the  observation 
equation  (50)  at  with  LiMobJ,  and  ACr  rto6^  set  to  zero  and  subtracting 

the  result  from  the  observed  value.  Next,  the  partial  derivatives  required  for  the  A 
matrix  are  defined.  Whitening  and  decorrelauon  of  the  observations  are  discussed  in 
the  FILTER  ALGORITHM  section. 
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Range  Partial  Derivatives 

As  mentioned  in  the  ESTIMATION  CONCEPTS  section,  the  partial  deriva¬ 
tives  of  the  observation  with  respect  to  the  parameters  are  required  for  the  Filter 
algorithm.  These  derivatives  are  obtained  by  differentiating  the  observation  equa¬ 
tion  with  respect  to  all  parameters  explicitly  present,  and  u»ng  the  chain  rule  and 
models  to  relate  them  to  the  solution  parameter  states.  Let  Pi>  =  ti  -r*,  pi,*  =  Ipi^l, 
and  R  =  Rijk,  then  the  partial  derivatives  of  range  with  respect  to  the  parameters  in 
the  observation  equation  are: 


3ri  8r*  P,,* 

dR  _  _  c 

aAi7  ”  ”  ~  Jo^ 

dR  _  1 

8AC|}^  sin  Ei^it 


(62) 

(63) 

(64) 


Using  these  partial  derivatives,  the  partial  derivatives  with  respect  to  the  solution 
parameters  are: 


dR 

dR 

dPj 

BKr, 

dR 

8R 

dr, 

ai\ 

dG, 

dR 

8R 

8ACfi/tj)  “ 

aACft/tp 

8R 

8R 

dA(, 

dCsv, 

8A(, 

acsv, 

dR 

dR 

dAt*  _ 

dC||^s^ 

8Ait 

acvs* 

dTi  2  dKRftj) 

(see  equations  (22)  and  (25)) 


ar,  2  dC,<tj) 

(see  equations  (29)  and  (32)) 


dR 

8AC«,  * 


dR  /  (hb*-toy 

8Aii  \  2 


to6»  “to  1 


dR  I 

(  to6»  “  to 

8aiA  \ 


1 


) 


(65) 

(66) 

(67) 

(68) 

(69) 


dR  _  dR  8r, 
8ei  8ri  de. 


(70) 
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dti 

where — is  interpolated  off  the  reference  trajectory  at 

d%i 


dR_ 


dR  _  ai 
8RP,  ~  dr,  dRP, 


flr* 

dASk 

dR 

=  —  MflCD)»Tk 

Ark 

(71) 

Tk  = 

(Ue  UN  Mv>* 

(72) 

uo 

II 

(73) 

Uv 

/  *«p  \  ^  “V 

=  [ykF  .  '“v  -  — 

(74) 

us 

UyXUo 
=  - 

luyXUol 

(75) 

Us 

=  UyXttc.  Un=  - 

UinI 

(76) 

3r. 

(77) 

ar 

where  — ^  is  interpolated  off  the  reference  tngectory  at  tok« 
BRP, 


dR  _  dR  dr, 
W,  ~  dTi  dT, 


(78) 


dTi 


0  if  toh»  s  tr, 

*/*r,  <  t,*.  s;  tr,  (79) 

[♦«/toAJ(^»7tr*.l^T/tr,)-^«j'tr,^^r/tra))]r  ifioh*  >  tr* 


where  4»«  =  partials  of  position  and  velocity  with  respect  to  epoch 

orbital  elements  obtained  by  interpolating  off  the  reference 
tri^ectory  at  the  appropriate  time 
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4>r,  =  partiaU  of  position  and  velocity  with  respect  to  the 

canonical  thrust  parameters  obtained  by  interpolating  off 
the  reference  tn^ectory  at  the  appropriate  time 

and  [  ]y.  denotes  the  first  three  rows  of  the  6X3  matrix 


SR 

dPM 


dTi  \dp  dq  d^t(T„>  aAt(To) )  dr*  dPdl 


(80) 


dTi  dr,  dr,  dr, 

where—.  — .  — r*  and -T-r— are  interpolated  off  the 

dp  dq  dA<fT„)  dAt(To) 

reference  trajectory  at  To  denotes  the  trajectory  epoch,  and 


ark 

dPM 


=  (BCD/  0 


0  —  Qftofc.— 

tok,  — 


**>•,  — y4>\ 


(81) 


dB  _  OR  dti 
d^  ”  ^  d^ 


(82) 


dr, 

where - is  interpolated  off  the  reference  trajectory  at 

dCC 


Observation  equations  and  partial  derivatives  for  all  other  data  types  are 
derived  from  the  range  observation  equation  and  partial  derivatives. 


RANGE  DIFFERENCE 

Range  difference  observations  can  be  either  of  two  types:  the  result  of  differ¬ 
encing  rann  observations  at  two  different  times  or  the  result  of  integrating  the 
D^pter-shifled  frequeucv  for  a  given  time  interval,  At/tp.  The  latter  type  includes 
differencing  accumulated  delta  ranges  that  are  continuous-count,  integ^ted 
Doppler  from  some  epoch.  Both  types  are  treated  the  same  in  the  measurement 
processing.  These  olmrvations  are  pairwise  correlated  if  the  end  of  one  range  differ¬ 
ence  interval  is  the  same  as  the  beginning  of  the  next  intervai.  This  correlation  is 
accounted  for  in  the  processing  as  describiM  in  the  FILTER  ALGORITHM  section. 
The  nonlinear  range  difference  (RD)  observation  equation  is  given  by 
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c  J 
10* 

( +  Ai,  —  (  a/,  (tuh»  —  At«/>>  +  At,(tu6,  —  At/ti>))j 

1 

i 

(83) 

— 

JO* 

\(  Ai*ftufc,)  +  ikikduht^)  —  ( Ai»  (tu*,  —  AtftoJ  +  Ai*(t„fc,  -  Atf(o)) 

1  • 

sin  E,,it(tuhs^  sin  Rt,t(tufut  “  At|j|>) 

where  tutu,  -  end  time  of  the  range  difference  interval 
and  Atfto  -  range  difference  interval  in  seconds. 

This  equation  is  just  the  range  observation  equation  evaluated  at  two  times  and 
differenced.  The  computed  value  is  the  difference  of  two  computed  ranm  values. 
Therefore,  the  partial  derivatives  are  the  range  partial  derivatives  differenced, 
where  the  derivatives  at  both  times  are  baaed  on  the  same  time  tj.  As  a  result  of  this 
differencing,  only  the  changes  in  the  clocks  over  the  interval  Atno  are  relevant,  i.e., 
dRD  _  BRD  _  ^ 

dAi,  dAii 

Differenced  and  doubly-differenced  range  dilTerence  data  types  are  not 
included  in  the  MSF/S  system.  However,  rai^  difference  data  can  be  processed  in  a 
mode  that  emulates  these  data  typos.  The  difierenced  range  difference  emulation  is 
based  on  the  idea  that  single  differencing  to  remove  the  satellite  clo<ik  fre^ency  off¬ 
sets  from  the  data  is  equivalent  to  solving  for  an  indeMndent  frequency  offset  for 
every  group  of  simultaneous  observations  of  that  satellite.  This  can  be  emulated  by 
setting  the  clock  frequency  offset  state  process  noise  variance  to  a  large  value  to 
decorrelate  the  estimates  between  mini-batch  steps.  The  satellite  frequency  drift 
state  should  not  be  solved  for.  For  doubly  differencing,  this  variance  adjustment 
must  also  be  done  for  all  station  clocks  except  the  master  cluck.  This  processing 
technique  fully  accounts  for  the  measurement  correlations  introduced  by  differ¬ 
encing.  However,  this  emulation  is  not  exact  if  more  than  one  observation  from  a 
given  satellite-station  pair  is  present  in  the  mini-batch  interval. 


DIFFEKKNCKU  RANGE 

Differenced  range  (DR)  observations  are  derived  by  differencii^  two  simulta¬ 
neous  observations  from  the  same  satellite  for  any  j>air  of  stations.  The  purpose  of 
this  differencing  is  to  eliminate  the  satellite  clock  from  the  observations.  No  account 
of  the  correlations  introduced  when  two  pairs  of  stations  have  one  station's  data  in 
common  is  included  in  processing  of  these  observations.  The  nonlinear  differenced 
range  observation  equation  is  given  by 
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DRi/Jk,4';  “ 

+  +  ^‘k(tob»))  —  (^tk‘  (tobs)  +  (84) 

^Cn  (t^bs)  ^^Ri^/^obs) 

+ - ! - ^ - 

sin  E,.*  ( t„bJ  sin  Ei^tf  ( tabs) 


where  k  and  k  are  indicies  specifying  the  two  stations.  Again,  the  computed  value  is 
obtained  by  di^erencing  two  range  computed  values,  and  the  partial  derivatives  are 
obtained  by  differencing  the  two  range  partial  derivatives.  Zeroes  are  used  for  the 
rang^e  partial  derivatives  involving  ^e  other  station.  This  differencing  results  in  the 
paruals  for  satellite  clock  states  being  zero. 

Range  data  can  be  processed  in  a  mode  that  emulates  this  data  type.  This 
emulation  is  based  on  the  idea  that  single  differencing  to  remove  the  satellite  clock 
time  offset  from  the  data  is  eouivalent  to  solving  for  an  independent  satellite  time 
offset  for  every  group  of  simultaneous  observations  of  that  satellite.  This  can  be  emu¬ 
lated  by  setting  the  clock  time  offset  state  process  noise  variance  to  a  large  value  to 
decorrelate  the  estimates  between  mini-batch  steps.  The  satellite  frequency  offset  and 
drift  states  should  not  be  solved  for.  This  processing  technioue  fully  accounts  for  the 
correlations  introduced  by  differencing.  However,  this  emulation  is  not  exact  if  more 
than  one  observation  from  a  given  satellite-station  pair  is  present  in  the  mini-batch 
interval. 


DOUBLY-DIFFERENCED  RANGE 

Doubly-differenced  range  (DDR)  observations  are  derived  by  differencing  two 
simultaneous  differenced  range  observations  from  any  pair  of  satellites  for  the  same 
pair  of  stations.  The  purpose  of  this  differencing  is  to  eliminate  the  station  clocks  from 
the  observations  in  addition  to  the  satellite  clocks  that  were  eliminated  by  the  first 
differencing.  No  account  of  the  correlations  introduced  by  this  differencing  technique 
is  included  in  the  processing  of  these  observations.  The  nonlinear  doubly-differenced 
range  observation  equation  is  given  by 


DDR(i^ii)jiif,it')  —  — rt(to5j()|  —  |r/t(,6s^  “r^i 

“  (|ri»rto6)i^  '^^k(^obs)\  ~\^i'(^obs) 

I 


if^obs)  f  ; 


—  ACft 


^•(tobs)  (  - - 


\sin  (^obtt)  bin  (^obtt) 

1 


\sin  Ei^b*  (^oba^  sin  Eit^b’  (^ubs^^ 


(85) 
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where  i,  i'  are  indices  specifying  the  two  satellites  and  k,  k  are  the  station  indices.  The 
computed  value  is  obtained  by  diiTerencing  two  differenced  range  computed  values, 
and  the  partial  derivatives  are  obtained  by  differencing  the  two  differenced  range 
partial  derivative  sets  with  appropriate  zeroes  for  irrelevant  parameters.  This 
differencing  results  in  all  clock  partials  being  zero. 

Range  data  can  be  processed  in  a  mode  that  emulates  this  data  type.  In  addition 
to  emulating  differenced  range  data,  all  station  clock  time  offset  states  (except  the 
master  clock)  should  have  their  process  noise  variances  set  to  large  values.  Also,  all 
station  frequency  offset  states  should  not  be  solved  for.  The  comments  given  under  the 
discussion  of  emulating  differenced  range  also  apply  to  doubly-differenced  range 
emulation. 


FILTKR/SMOOTHER  PROCESSING  FLOW 


Figure  4  gives  a  functional  definition  of  the  processing  flow  within  the  MSF/S 
system  of  programs.  The  Filter  is  initialized  at  to  or  restarted  at  an  arbitrary  mini¬ 
batch  step  if  previous  filtering  has  been  done.  At  each  mini-batch  step  a  measurement 
update  is  performed  first,  followed  by  solving  the  equations  and  generating  diagnostics 
if  required,  and  then  propagating  to  the  next  mini-batch  step.  No  propagated  solution 
is  ever  computed  unless  observations  are  not  present  in  a  given  mini-batch  interval 
and  solutions  at  each  mini-batch  step  are  required.  This  process  is  repeated  until  the 
solution  and  diagnostics  at  t\  are  completed.  At  this  point  the  y  parameter  solutions 
are  final,  i.e.,  no  smoothing  of  the  y  parameters  is  possible.  If  stochastic  orbit-related 
parameters  are  not  present,  the  resulting  orbital  element  and  constant  force  model 
parameter  corrections  can  be  applied  to  their  initial  values  and  used  to  reinte^ate  a 
trajectory  if  the  orbit  is  not  converged,  or  an  improved  trajectory  can  be  linearly  propa¬ 
gated  as  in  a  batch  fit.  The  polar  motion  and  station  coordinate  tables  are  updated  at 
&is  point.  The  updated  station  coordinates  are  required  for  the  residual  generation 

Erocedures.  The  appropriate  information  arrays  and  partial  derivative  matrices  must 
e  saved  from  the  Filter  at  each  mini-batch  step  to  be  used  in  the  smoothing  process.  A 
Filter  propagated  trajectory  can  be  created  at  tnis  time  if  required  and  also  a  set  of 
improved  initial  conditions. 

Two  paths  within  the  Smoother  are  possible.  If  state  and  covariance  estimates 
are  required  the  left  side  is  followed.  If  state  estimates  only  are  required  the  right  side 
is  followed.  In  the  first  case  the  smoothing  arrays  are  initialized  at  t^  and  smoothing 
proceeds  in  reverse  time  order.  At  each  mini-batch  step  smoothing  arraj^s  are 
manipulated  followed  by  solving  the  equations  and  generating  diagnostics.  This  is 
repeated  at  each  step  until  the  process  terminates  at  to.  Smoothed  trajectories  can 
then  be  propagated  and/or  the  SATRACK  covariance  matrices  can  be  computed.  In  the 
second  case  the  state  estimates  are  initialized  at  ts  as  in  the  first  case  and  smoothing 
proceeds  in  reverse  time  order.  At  each  mini-batch  step  the  state  estimate  is  computed 
based  on  the  estimate  at  the  previous  step,  the  information  arrays  saved  in  the  Filter, 
and  the  stochastic  state  equations.  This  process  is  repeated  until  to  is  reached.  No 
covariance  information  can  be  computed  for  this  option  but  the  rest  of  the  diagnostics 
can.  Trajectories  can  then  be  propagated  followed  by  generation  of  residuals. 


36 


NSWCTR  87-187 


FIGURE  4.  FILTER/SMOOTHER  PROCESSING  FLOW 
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The  rest  of  this  report  will  detail  the  square  root  information  implementation 
of  the  Filter  and  Smootlier  algoriUims,  the  method  of  obtaining  the  solution  and 
diagnostics  (these  are  identical  in  each  algorithm  except  for  handling  the  y  parame¬ 
ters  and  for  residual  generation),  the  trajectory  propagation  procedures,  and  the 
SATRACK  covariance  matrix  generation  formulation.  Derivations  of  some  of  the 
equations  and  showing  that  they  are  equivalent  to  the  standard  Kalman  (Ilter^TS 
smoother  equations  are  included  in  the  appendicies.  Also  a  discussion  of  properties 
of  Householder  transformations  and  the  upper  triangular  matrix  inversion  method 
are  given  in  the  appendicies, 


FILTER  ALGORITHM 

The  Filter  can  be  described  as  a  square  root  information  Hlter  based  on  the 
form  of  the  state  equations  given  in  the  STATE  EQUATIONS  section  and  the 
observation  equations ^ven  in  the  OBSERVATION  EQUATIONS  AND 
PARTIAL  DEVIRATI VES  section.  The  Filter  consists  of  four  processing  steps, 
the  last  three  of  which  are  repeated:  Initialization/Restart,  Measurement  Update, 
Solution  and  Diagnostics  (optional) ,  and  Propagation.  After  initialization  or  restart 
of  the  information  arrays,  at  each  mini-batch  step ,  tj,  a  measurement  update  is  per¬ 
formed  if  there  are  any  observations  in  the  corresponding  mini-batch  interval.  Then 
the  solution  and  diagnostics  are  computed  if  required  ana  the  information  array  is 
augmented  and  propagated  from  tjtatj+f  Portions  of  the  information  array  required 
for  smoothing  are  then  saved.  This  cycle  repeats  until  a  measurement  update  at  tyv 
has  been  done.  The  y  parameter  information  at  this  time  is  final  and  is  used 
throughout  the  smoothing  procedure.  Details  for  each  processing  step  are  given 
below. 


INITIALIZATION/RESTART 

The  primary  purpose  of  this  step  is  to  set  up  and  initialize  the  information 
arrays.  The  parameters,  if  present,  are  always  ordered  as  given  in  the  state  equation 
description,  i.e.,  p  parameters  with  orbit-related  ones  before  measurement-related 
ones,  X  parameters,  and  then  y  parameters.  In  addition  since  both  orbit-related  p 
parameters  and  tropospheric  refraction  are  modeled  as  first-order  Gauss-Markov 
processes,  these  are  at  the  top  of  the  parameter  list  so  that  the  propagation  step  com¬ 
putations  can  be  made  more  efficient.  The  following  notation  is  used  in  the  rest  of 
this  report. 

No  =  number  of  orbit-related  p  parameters,  maximum  of  6Nsv 

Nm  =  number  of  measurement-related  p  parameters,  maximum  of 
3Nsv  +  3Nms 

Np  =  number  of  p  parameters  ( must  always  be  >  0)  =  No  +  Nu 
Nom  =  number  of  Gauss-Markov  p  parameters,  maximum  of  6Nsv  + 

N,  =  number  of  x  parameters,  maximum  of  6Nsv 

Ny^  =  number  of  station- related  y  parameters,  maximum  of  3iVj|S 
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Sy  =  number  of  orbit-related  y  parameters,  maximum  of 
3Nsv  +  3Nr  +  Npm  + 

Ny  =  number  of  y  parameters  =  N,  +  Ny^ 

SroT  =  total  number  of  parameters  =  Np  +  Nx  +  Ny 
where  Nsv  =  number  of  satellites 
Njus  =  number  of  stations 
Nr  =  number  of  thrusts 
NpM  =  number  of  polar  motion  parameters  =  3 
Ngc  —  number  of  gravity  coefficient  parameters 

The  parameter  set  for  a  particular  fft  is  selectable  but  with  the  following 
restrictions: 

1 .  All  orbit-related  parameters  except  for  thrusts,  polar  motion,  and  gravity 
coefficients  are  present  for  all  satellites  and  individual  parameters  must  be  selec¬ 
tively  deweighted  if  required. 

2.  All  measurement-related  parameters  are  present  for  all  satellites  or  sta¬ 
tions  and  individual  parameters  must  be  selectively  deweighted  if  required.  The 
clock  parameters  for  the  master  station  are  automatically  deweighted  and  the 
corresponding  white  noise  spectral  densities  are  set  to  approximately  zero. 

3 .  All  station-related  parameters  are  present  for  all  stations  and  must  be 
selectively  deweighted  if  required.  If  orbits  and  station  coordinates  are  being  solved 
for  simultaneously,  the  east  component  of  one  station  should  be  deweighted. 

4.  If  only  differenced  range  and  doubly-differenced  range  are  being  processed, 
satellite  clock  parameters  are  not  included  in  the  state  equations.  If  only  doubly- 
differenced  data  are  being  processed,  satellite  and  station  clock  parameters  are  not 
included. 

5.  If  only  range  difference  data  is  being  processed,  all  satellite  and  station 
time  offset  parameters  should  not  be  deweighted  unless  all  the  white  noise  spectral 
densities  are  zero. 

6.  Stochastic  orbit-related  parameters  can  only  be  solved  for  if  pseudoepoch 
orbital  elements  are  also  present  in  the  equations. 

The  following  notation  is  used  in  the  rest  of  this  report: 

'  indicates  a  Filter  predicted  quantity 

*  indicates  a  Filter  estimated  quantity 

*  indicates  a  Smoother  estimated  quantity 
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The  initial  forms  of  the  two  information  arrays  at  to  are  given  by 


p-x  information  array 


Ntot  1 


y  information  array 

iK  »)  J  N, 

Ny  1 


(86) 


where  is  an  NpX  Np  diagonal  matrix  with  each  element  of  the  form  1/Op  and  is 

an  Ny  X  Ny  diagonal  matrix  with  each  element  of  the  form  1/Oy.  and  Oy  are  a  priori 

parameter  sigmas  in  internal  units.  (The  y  array  is  not  present  if  Ny  =  0.)  The 
matrix  is  N^XN,  and  has  the  form 


where  each  is  a  6X6  matrix  computed  for  each  satellite  by 


ft.  =  CkacR^T' 


(88) 


/' 


2 

Or 


where  Crac  = 


\ 


2 

Oa 


0 


Oc 


Or 


0 


OA 


oc 


V 


(89) 


o's  are  inj^ut  a  priori  si^as  on  radial,  along-track,  and  cross-track  position  (in  km) 
and  velocity  (in  km/sec)  at  the  fit  epoch. 


R  0 


R  = 


(90) 


^  .  r  X  V  .  r  X  V 

R  =  r  ^  X  r 


(91) 
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r,  V  are  inertial  coordinates  of  the  satellite  at  the  Ht  epoch  interpolated  off  of  the 
reference  trajectory. 


/ 

dx 

dx 

1  da 

de^ 

dCl 

dy 

<iy 

dy 

r  = 

da 

dfit 

dCi 

dz 

dz 

diz 

\  da 

deg 

dCl 

\ 


(92) 


T  is  a  matrix  of  partial  derivatives  of  position  and  velocity  at  the  fit  epoch  with 
respect  to  orbital  elements  at  the  trajectory  epoch  interpolated  off  of  the  reference 
trajectory. 


Rt  is  just  the  fit  epoch  RAC  position  and  velocity  sigmas  transformed  to  a  covariance 
matrix  on  orbital  elements  at  the  trajectory  epoch,  inverted  and  with  the  square  root 
taken. 


To  restart  the  Filter  the  assumption  is  made  that  all  physical  constants, 
spectral  densities,  decorrelation  times,  and  other  quantities  must  be  the  same  as 
those  used  in  the  last  execution  of  the  Filter.  The  two  information  arrays  are  then 
initialized  with  the  results  of  the  propagation  step  from  tt-i  to  tf  as  follows: 


p-x  information  array 


y  information  array 


ftpx 


(Ry 


(93) 


MEASUREMENT  UPDATE 

All  observations  in  the  mini-batch  interval  (^tj  —^•*7+  are  processed 

simultaneously  because  computationally  this  is  the  most  efficient.  Measurement 
updating  is  carried  out  by  augmenting  the  propagated  information  arrays  with  the 
whitened  and  deccrrelated  measurement  matrix  and  residuals  (A  z)  and  tranforming 
this  expanded  array  using  Householder  orthogonal  transformations.  The  details  of 
these  procedures  are  given  as  follows: 

At  At 

1.  All  observations  with  observation  times,  t„6»,  such  that  tj  -  tj  +  — 

are  processed  for  the  stations  and  satellites  selected  unless:  2  ^ 

a.  this  data  type  is  not  being  processed  or  t„fr«  lies  outside  of  the  subspan 
for  this  particular  data  type, 

b.  the  corresponding  pass  number  is  in  the  list  of  passes  to  be  deleted. 
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c.  the  observation  sigma  indicates  that  this  observation  has  been 
previously  tagged  (negative  sigmas  indicate  this), 

d.  anv  elevation  angle  is  less  than  B^a.  -cm  input  tolerance  in  de^ees,  or 

e.  a  clock  event  for  we  given  station  or  satellite  is  present  in  this 
interval. 

At  the  start  time,  the  reduced  mini-batch  step  span  the  mini-batch 

/  At 

interval  is  defined  as  Us'***'  -  +  ^  |.  At  the  end  time,  of  the 

'  it  • 

psp 

]. 


tg 


Red. . 


reduced  mini-batch  step  span  the  mini-batch  interval  is  defined  as 

^iRed.  At^ 

If  no  observations  are  present  in  a  given  mini-batch 


interval,  the  Householder  transformation  is  still  done  in  order  to  upper  triangularize 
the  array,  i.e.,  the  transformation  indicated  in  equation  (103)  below  is  done  with  the 
(A  z)  rows  absent.  If  no  solution  is  required  the  processing  skips  to  the  propagation 
step. 


2.  For  all  remaining  observations,  the  partial  derivative  matrices 
(designated  A)  relating  the  observations  to  each  parameter  are  computed  along  with 
the  residuals  (designated  O-C')  as  described  in  the  OBSERVATION  EQUATIONS 
AND  PARTIAL  DERIVATIVES  section.  Only  those  partial  derivatives  with 
respect  to  the  relevant  satellite  and  station  parameters  are  non-zero.  All  others  are 
zero  with  each  column  corresponding  to  a  parameter  as  ordered  in  the  state 
equations.  These  quantities  are  denoted  by 

A=fA'pA,  A*,)  z  ^OC 

Mj  X  (Np  +  N,  +  Ny)  Mj  X  J 

where  Mj  =  number  of  observations  in  the  mini-batch  interval  centered  at  tj. 

3  All  obsei  vations  are  then  whitened  and  decorrelated.  To  whiten  ran^^e, 
differenced  range,  and  doubly-differenced  range  each  row  of  the  partial  derivative 
matrix  A  and  the  residuals  z'  is  divided  by  its  observation  sigma,  o„6ti>  or  by  Oubc^w 
that  data  type  if  Oo6ii<  o„6,  to  get  A  and  x.  For  the  range  difference  data  where 
tobs  —  *  toks  (^  viou^,  i.e.,  the  beginning  of  the  ran^e  difference  interval  is  not 

identical  to  the  rvation  time  of  the  previous  observation  from  this  same  station- 

satellite  pair,  a  similar  procedure  is  used.  However,  if  the  times  do  coincide  the 
observations  are  correlated  and  must  be  whitened  and  decorrelated  before  being 
processed.  The  range  difference  data  are  whitened  and  decorrelated  as  follows  (see 
Appendix  F  for  the  derivation  of  this  procedure): 

If  the  observation  is  not  tag^d  the  two  input  sigmas  for  the  observations 
differenced  to  obtain  the  range  difference  observation,  designated  oahr„.,  nnd  oadr^, 
are  possibly  redefined  as  follows: 


vf 

OHOm/n 

A  V  <  OadK„  ,  <  — ■== 

V2 

(95) 

oadr^ 

VT 

11  <  - ■= 

VT 

(96) 
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The  range  difference  observation  sigma  is  then  given  by 


Ord„ 


^ADH„ 


(97) 


If  oadr  ,  =  0  or  the  last  range  difference  observation  for  this  satellite-station  pair 
was  nol  processed,  the  observation  is  whitened  by  dividing  Ard^  and  zhu  by  ord  . 
5„  is  then  set  equal  to  oro,,,  and  5,„  Aru  ,  and  xru^  are  saved.  If  Oao«„  ,  *  6  tne 
observation  is  whitened  and  decorrelated  as  follows; 


_  -  OaUR 

Pn-i 

Snl 

(98) 

On 

=  (Or/)„- 

pU)* 

(99) 

Aro„ 

=  (Ard„  - 

Pn  l  ^RD^  JI^n 

(100) 

ZRD„ 

=  (zrd„  - 

Pn-1  Zro^  ^)IOn 

(101) 

On,  Aro^,  and  zrq^  are  then  saved.  It  is  assumed  that  if  an  entire  mini-batch  interval 
is  processed  before  the  same  satellite-station  pair  occurs  again,  the  new  observation 
is  uncorrelated  with  the  previous  observation.  This  is  done  so  that  information  does 
not  have  to  be  saved  for  more  than  one  mini-batch  interval. 


The  effect  of  the  accumulated  clock  noise  during  the  interval  is  not 

modeled  so  the  minimum  observation  sigma  for  range  difference  data  should  be  used 
to  account  for  this  effect. 


4.  The  propagated  p-x  and  y  information  arrays  without  the  terms  required 
for  smoothing  are  given  by 


(Np  +  N,)  X  (Np  +  N^  +  Ny -hi) 


(Ry  Zy)j  —  (Ry  Zy)j  j 


(102) 


Ny  X  (Ny  +  1) 


The  measurement  update  is  done  in  two  steps  if  y  parcuneters  are  present.  The 
p-x  information  array  is  augmjented  by  (A  z),  ana  a  sequence  of  Householder 
orthogonal  transformations,  fp^y  are  applied  to  zero  out  elements  below  the  diagonal 
of  the  first  Np+Nj  columns: 
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rNp  +  N,  +  Alj>  X  rNp  +  N,  +  Ny+i; 


(103) 


where  Rp  and  R,  are  upper  trian^lar  matrices.  If  no  y  parameters  are  present  the 
columns  with  a  y  subscript  are  absent,  z  =  «,  and  the  measurement  update  is 
complete.  The  e  term  is  related  to  the  sum-of-squares  of  weighted  residuals.  If  y 
parameters  are  present,  t)ie  measurement  update  is  completed  by  augmenting  the  y 
information  arrav  with  {Ay  i)  and  applying  anoUier  sequence  of  Householder 
orthogonal  transformations,  f  to  zero  out  elements  below  the  diagonal  of  the  first 
Ny  columns: 


A 

T 


y 


(104) 


(Ny+Mj)  X  (Ny-hl) 

where  both  Ry  and  Ry  are  u^er  trianmilar  matrices,  i.e.,  the  y  information  array  is 
always  upper  triangular.  (Ry  ?^)  could  be  added  as  extra  rows  in  eouation  (103)  above 
but  this  is  not  done  to  save  storing  a  large  block  of  zeroes  that  would  never  change. 

The  transformation  matrices  f  p,  and  fy  are  not  explicitly  computed. 
Householder  transformations  can  be  carried  out  by  operating  on  the  columns  of  a 
matrix  one  at  a  time  as  follows: 


Let  R  be  an  arbitrary  m  Xn  information  array.  To  zero  out  all  elements  below 
the  diagonal  of  the  Hrst  column  define  a  scalar  s  and  a  vector  u  by 


s  = 


/  ^ 

-sgn(R(l,l))f  S  [Rfi.Dp' 


(105) 


u(l)  =  R(I, l)-s 
u(i)  =  R(i,l)  i  = 

1 

a  =  - 

3U(1) 


(106) 

(107) 

(108) 
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m 

Then  define  Yj  =  a  S  u(0  R({  J)  (109) 

f  =  J 


The  effect  of  the  transformation  on  a  column  /  is  then  written  as 

T„(R(i,/))  =  R(i,i)  +  YjU(l)  i  =  l,2...,m  (110) 


This  results  in  the  first  column  being 


as  shown  in  Appendix  G.  This  is  the 


intended  result  and  this  must  be  done  to  every  column  of  the  array.  This  procedure 
is  then  repeated  to  zero  out  Gie  below  diagonal  elements  of  the  second  column  with  u 
deHned  by  the  last  m  —  1  rows  of  this  column.  This  procedure  (in  which  the  matrix 
operated  on  decreases  in  both  dimensions  at  each  step)  is  repeated  until  the  required 
number  of  columns  have  been  transformed.  A  property  of  Householder  transforma¬ 
tions  t^at  affected  the  implementation  of  equations  (103)  and  (104)  is  discussed  in 
Appendix  G  also. 


5.  If  any  diagnostics  are  to  be  computed  or  a  trmectory  is  to  be  propagated,  the 
solution  is  then  computed  as  described  in  the  SOLUTION  AND  DIAGNOSTICS 
section.  If  a  maqter  station  switch  has  taken  place  in  this  interval  the  two  diagonal 
elements  of  tihe  Rp  array  corresponding  to  the  new  master  station  clock  parameters 
are  set  to  large  values  after  solution.  If  the  last  mini-batch  interval  centered  at  ts 
has  been  processed  the  filtering  is  complete  and  a  solution  is  computed  for  the  y 
parameters  if  present. 


PROPAGATION 


Propagation  from  t,  to  involves  modifying  the  p-x  information  array  to 

incorporate  the  effeqts  of  process  noise.  Bias  parameters  are  unaffected  by  process 
noise  so  that  =  Propagation  also  involves  generation  of  auxiliary  arrays 

required  for  smoothing.  The  propagation  step  is  based  on  the  state  transition  and  Rw 
(derived  from  Q)  matrices  described  in  the  STATE  EQUATIONS  section.  An 
augmented  p-x  information  array  is  upper  triangularized  over  the  first  Np  columns  to 
carry  out  the  propagation  as  follows: 


(111) 
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where  R»  is  an  upper  triangular  matrix  and  *  denotes  the  smoother-related  matrices 
required  to  obtain  smoothed  state  and  covariance  estimates  at  tj.  The  Househoulder 
transformations  are  done  as  described  above  in  the  Measurement  Update  processing 
step  except  that  a  computational  reduction  is  made.  All  columns  of  the  —  matrix 

have  zeroes  below  the  diagonals  down  to  and  includina  row  Np.  When  applying  the 
transformation  to  the  first  column,  rows  2  through  Np  for  all  columns  do  not  change 
because  of  a  property  of  Householder  transformations  given  in  Appendix  G.  These 
rows  are  ignored  in  the  computations  to  save  on  computer  time.  This  reduction 
applies  to  the  first  N.  - 1  columns  as  they  are  zeroed  out  below  the  diagonal.  Another 
property  of  Householder  transformations  given  in  Appendix  G  is  also  applied  in 
implementing  equation  (111)  and  is  possible  because  Ru,  is  diagonal  for  tne  first 
No'm  rows  and  columns.  When  processing  the  first  column  the  only  additional  column 
that  changes  out  uf  the  the  first  2  Np  columns  is  the  Np+1  column.  When  processing 
the  second  column  only  the  Np+1  and  Np+2  columns  change  out  of  the  first  2  Np 
columns.  One  additional  column  is  affected  until  column  Ngm+  1  is  processed.  The 
part  of  Rp  corresponding  to  the  Gauss-Markov  parameters  is  upper  triangular  and 
the  part  of  R^  corre^nding  to  these  same  parameters  is  lower  triangular.  This  is 
the  reason  why  the  Uauss-Markov  parameters  are  ordered  first  in  the  state  equa¬ 
tions.  The  Ru,  and  M  arrays  are  a  function  of  the  mini-batch  step  size  so  they  are 
adjusted  starting  at  and  then  reset  to  their  original  values  at  .  Portions  of 

the  Ru,  matrices  corresponding  to  the  satellite  and  station  clocks  are  adjusted 
appropriately  for  events  as  reouired.  The  predicted  solution  is  never  computed  in  a 
square  root  information  Hlter  because  either  a  full  matrix  inversion  or  a  House¬ 
holder  transformation  followed  by  inversion  of  an  upper  triangular  matrix  is 
required. 

Appendix  H  discusses  the  derivation  of  the  measurement  update  equations 
(103)  and  (104)  and  the  propagation  equations  (111).  Also  this  appendix  shows  the 
mathemaUcal  equivalence  between  Uiese  equations  (actually  a  more  general  form  of 
the  propagation  equations)  and  the  Kalman  filter  equations  given  by  equations 
(3)  tJiru  (7)  in  the  ESTIMATION  CONCEPTS  section. 


SMOOTHER  ALGORITHMS 


Two  smoother  algorithms  are  included  The  first  uses  Householder  ^ans- 
formations  to  upper  triangularize  smoothing  information  arrays  from  which  both 
state  and  covariance  estimates  can  be  computed.  This  is  referred  to  as  the  Array 
Smoother  here  but  is  commonly  called  the  square  root  information  smcMther.  The 
second  combines  the  stochastic  state  equations  and  data  equations  derived  from  the 
information  arrays  saved  in  the  propagation  step  of  the  Filter  to  obtain  state  esti¬ 
mates  only,  i.e.,  no  covariance  information  is  available.  This  is  referred  to  as  the 
State  Only  Smoother.  The.  Array  Smoother  consists  of  three  processii^  steps,  the 
last  two  of  which  are  repeated:  Initialization,  Array  Smoothing,  and  Elution  and 
Diagnostics.  At  each  mini-batch  step  t,  a  smoothing  array  is  constructed  using  the 
array  solved  at  /  and  the  array  saved  in  the  Filter  when  propagating  from  tj  to 
tj  +  i.  This  array  is  upper  triangularized  for  2Np -f  N,  columns  and  combined  with  the 
y  parameter  solution  matrices  from  (which  remain  constant  throughout  the  span) 

to  obtain  the  state  and  covariance  estimates.  This  cycle  repeats  until  the  solution  at 
to  has  been  computed.  Details  for  the  first  two  processing  steps  are  given  below 
followed  by  a  description  of  the  State  Only  Smoother. 
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INITIALIZATION 


1‘he  purpose  of  this  step  is  to  set  up  and  initialise  the  information  arrays.  This 
step  is  common  to  the  two  smoother  almrithms.  The  order  of  the  parameters  is  the 
same  as  that  used  in  the  Filter.  At  ts  tne  arrays  are  initialised  with  the  arrays 
determined  after  the  measurement  update  at  done  in  the  Filter.  The  Smoother 
solution  at  ts  is  identical  to  the  Filter  solution  at  ts-  The  information  arrays  at  ts 
are  given  by 

p-x  information  array  y  information  array 


Np  N,  Ny  1 


Ny  1 

(Ry  iy)/^  }  Ny 


(112) 


where  R^,  R„  and  R.  are  upper  triangular  matrices.  The  y  information  array  does  not 
change  m  the  smoothing  process. 


ARRAY  SMOOTHING 


Smoothing  of  the  p-x  information  array  at  time  L  is  accomplished  by  applying  a 
sequence  of  Householder  transformations  to  sero  out  elements  below  the  diagonal  of 
the  first  2Np + N,  columns  of  an  augmented  p-x  information  array  as  follows: 


T 


Np  N,  Ny 

Rp(t>)  +  Spp(tj)M  +  R^(ty)Vp^  Rp,(tj) 

R*p{tj^,)M  +  R^(tj^t)Vp^  R^itj^t) 

R*(tj+/)Vpj  RJftj+j)  RJyCt^+i) 


J 


(113) 


/Rpp(t,)  R'pitj)  R'p,(tj)  R'pyitji  g'pdji 
0  Rp(tj)  Rpx(t^)  Rpy{tj)  Xp{tj) 

0  0  Rl(tj)  Rlyitj)  mtj) 


where  the-  terms  correspond  to  the  *  terms  saved  in  propagating  in  the  Filter  from  tj 
to  tj+i.  Rp,  Rpp,  RJ^,  and  r;  are  all  upper  triangular  matrices.  The  Householder  trans¬ 
formation  is  earned  out  as  described  above  in  the  Measurement  Update  processing 
step.  These  computations  also  take  advantage  of  the  sparseness  or  the  Vp^  matrix  and 
the  fact  that  M  is  diamnal.  The  terms  identified  with  the '  superscript  are  not 
required  for  any  further  computations.  The  solution  method  for  Uiis  approach  is 

fiven  in  the  SOLUTION  AND  DIAGNOSTICS  section.  Appendix  I  contains  a 
erivation  of  these  smoothing  equations  and  also  shows  the  mathematical  equiva¬ 
lence  between  these  equations  (actually  a  more  general  form)  and  the  RTO  smoother 
equations  given  by  equations  (8)  thru  (10)  in  the  ESTIMATION  CONCEPTS 
section. 
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STATE  ONLY  SMOOTHING 

This  approach  utilises  the  stochastic  state  eauations  involving  the  p  and  x 
parameters  along  with  the  data  equations  derivea  from  the  information  arrays  saved 
in  the  Filter  propantion  step  to  recursively  generate  smoothed  estimates  for  the  p 
and  X  parameters.  The  initialisation  arrays  given  above  are  solved  at  u  to  get  Apjy, 
Axjv.  and  Ay^.  Then  given  the  smoothed  state  solution  at  r,  the  stnootned  state 
solution  at  tj  Is  just  given  hy 


Ap;  =  [ilp(tj)l  '  [Vtj)-ilpp(t^)Ap,%,-ftp,(tj)Ax;+i~ftp,(t;)Ay;i 
Ax;  =  Ax;+,-Vp^A|>’ 


(114) 

(116) 


where -denotes  matrices  saved  in  propagating  from  and  Rpitj)  is  always 

upper  triangular.  All  elements  of  V„  which  multiply  non-orbit-related  p  parameters 
are  sero  and  orbit-related  p  parameters  for  a  given  satellite  affect  the  x  parameter 
solution  for  that  satellite  onV>  This  procedure  requires  the  inverse  of  an  N.X  N. 
upper  triangular  matrix  at  each  mim-batch  step  instead  of  an  (Np-f  N,)X(Np+N,) 
matrix  as  in  the  Array  Smoother  algorithm. 


SOLUTION  AND  DIAGNOSTICS 


The  solution  method  and  the  diagnostic  computations  are  almost  identical 
between  the  Filter  and  Array  Smoother  and  are  therefore  described  here  together. 
The  state  solutions  and  all  diagnostics  depend  on  first  computing  the  inverses  of  the 
upper  triangular  p-x  and  y  information  arrays.  These  inverses  are  always  computed 
in  the  Array  Smoother  but  are  only  computed  in  the  Filter  if  any  diagnostics  are 
required  or  a  trajectory  is  to  be  propagated.  The  possible  diagnostics  are  correlation 
coefficients  (at  every  nth  mini-batch  step),  transformed  corrections  and  standard 
deviations,  total  clock  offsets,  and  residuals  and  signal-to-noise  after  fit.  For  ^e 
State  Only  Smoother  all  diagnostics  are  available  except  for  the  standard  deviations 
and  correlation  coefficients. 


SOLUTION 


The  solution  (state),  AX^,  and  covariance,  Pj,  estimates  are  computed  as  follows: 
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f  ^  T 

Rj  f<j 


/Rp  Rp: 

\0  R^ 


Rp  Rpx\  /Rpy 

0  R,Ar,v, 


(117) 


(118) 


where  all  quantities  denoted  *  in  the  Filter  and  *  in  the  Array  Smoother  are 
evaluated  at  tj  except  for  the  Smoother,  in  which  R^'and  Zy  are  fixed  at  their  values 
at  t/v.  All  Rp,  Rx,  and  R,  matrices  are  always  upper  triangular  so  that  both  ( Rp  Rp*\  * 

\o  rJ 

and  Ry  are  upper  triangular  also,  since  upper  traingularity  is  preserved  by  inversion. 
The  inversion  of  upper  triangular  matrices  is  discussed  in  Appendix  J.  The  solution 
is  always  computed  for  tlie  y  parameters  at  in  the  Filter  and  may  be  computed  at 
each  mini-batch  step  if  the  y  parameter  only  or  full  diagnostics  are  being  computed. 
The  solution  information  is  saved  in  the  Filter  only  if  a  trajectory  is  to  be  later 
propagated.  In  the  Array  Smoother  this  information  is  always  saved  because  it  is 
required  for  trajectory  propagation  and  residual  generation.  The  full  covariance 
matrix  is  computed  only  for  mini-batch  steps  for  which  the  correlation  coefficient 
matrix  is  required.  Certain  submatrices  oi  the  covariance  matrix  are  required  for 
deriving  the  standard  deviations  on  transformed  corrections. 

At  the  last  mini-batch  time  ts  in  the  Filter  the  final  y  parameter  solutions  are 
available.  Each  set  of  coordinates  in  the  station  coordinate  table  is  then  updated  as 
follows,  if  these  parameters  were  improved: 

/  \  /\  /  m  (l-e^sin^4>)"^AEN  \ 

/  \  /  \  /  n  QuankCOS^  \ 

/  n  aKarih(i-e^) 


Updated 


Updated 


ysF  I  +  T 


where 


/AE\ 

(AN 

\AV/ 


(120) 


are  the  coordinate  corrections  at  In  and  T  is  the  transformation  matrix 


defined  in  the  OBSERVATION  EQUATIONS  AND  PARTIAL  DERIVATIVES 
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section-  equations  (72)-(76).  Each  daily  entry  in  the  polar  motion  table  is  also 
updated  as  follows: 


where  t  —  to  is  in  seconds. 


(121) 


Improved  initial  conditions  at  the  trajectory  epoch  To  can  be  computed  in  the 
Filter  if  no  stochastic  orbit-related  parameters  are  present.  These  improved  initial 
conditions  would  primarily  be  used  to  reintegrate  a  trajectory  if  convergence  has  not 
occurred  when  emulating  a  batch  fit  processor.  The  required  equations  are  given  as 
follows: 


where  the  partial  derivatives  olFj>osition  and  velocity  with  respect  to  orbital  elements 
are  obtained  by  interpolating  oft  of  the  trajectory  at  Tq. 


-I-  AT, 


i  =  l,2 


(123) 

(124) 


COhRELATlON  COEFFICIENTS 


The  correlation  coefficient  matrix  is  computed  every  nth  mini-batch  step  in 
terms  of  the  actual  solution  parameters  and  not  a  transformed  set.  Each  element  of 
the  covariance  matrix  is  computed  as  follows: 


Ntot 

Pm,n  =  Pn.m  =  S  n,m  =  1,2,...,NtOT 

t  —  maxim/i) 

where  reR  '  and  peP  =  R  'R  ^ 

and  {=  max(m/i)  takes  advantage  of  the  upper  triangular  form  of  R  '.  Each 
correlation  coefficient  is  then  computed  as  follows: 


Cm,n  —  C/».m  — 


J. 

Pm.n 


^Pm,m 


if  m=n 
if  m^n 


Pn.li 


n,m  —  1,2.,.,Ntot 


(125) 


(126) 
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Only  a  lower  trian^lar  array  of  correlation  coeHlcients  is  computed  since  the 
covariance  matrix  is  symmetric.  The  correlations  between  any  pair  of  y  parameters 
in  the  Smoother  is  the  same  at  every  mini-hatch  step. 


TRANSFORMED  CORRECTIONS  AND  STANDARD  DEVIATIONS 


All  solution  states,  as  mentioned  before,  are  corrections  to  nominal  values  in 
internal  units.  On  option  these  corrections  and  their  corresponding  covariances  are 
converted  to  more  meaningful  corrections  and  standard  deviations  before  being 
printed  or  plotted.  If  the  state  only  smoothing  option  is  selected  no  standard  devia¬ 
tions  can  be  computed,  y  parameter  only  corrections  and  standard  deviations  can 
also  be  computed  in  the  Filter. 


Stochastic  radiation  pressure  parameter  and  all  orbit-related  y  parameter 
corrections  and  standard  deviations  are  unchanged.  Gravitational  acceleration 
parameter  corrections  and  standard  deviations  are  converted  to  10  '^km/sec^  for 
plotting  so  that  they  are  in  the  same  units  as  the  y-axis  acceleration  parameter. 
Tropospheric  refraction  parameter  corrections  and  standard  deviations  arc  con¬ 
verted  to  cm  for  plotting.  Satellite  clock  parameter  corrections  and  covariances  are 
converted  from  pseudoepoch  state  to  current  state  representations  as  follows: 


Aio  1  =  ^C'sy^^SvUj) 

V“'J 


(127) 


PCsv(^j)  -  ^Csv 


(128) 


Frequency  drift  terms  are  converted  from  ppm/sec  to  parts  in  10‘Vday,  frequency 
offset  terms  are  converted  from  ppm  to  parts  in  10’^  and  time  offset  terms  are  con¬ 
verted  from  psec  to  nsec.  Station  clock  parameter  corrections  and  covariances  are 
also  converted  from  pseudoepoch  state  to  current  state  as  follows: 


(129) 


PCys^^j)  - 


(130) 


Frequency  offset  terms  are  converted  from  ppm  to  parts  in  10'-  and  time  offset  terms 
are  converted  from  psec  to  nsec.  Only  the  square  roots  of  the  diagonals  of  Pcgy  and 
P(^Ms  required. 
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For  each  satellite  the  pseudoepoch  state  orbital  element,  radiation  pressure, 
thrusUs),  and  gravity  coefHcient  corrections  and  covariances  are  converted  to 
position  and  velocity  corrections  and  covariances  in  the  RAC  reference  frame  as 
follows: 


4‘RpARPj  +  4>t, At +  4>pmAPM,  +  ^gcAGC^  ) 


(131) 


PRAc(tj)  =  (^e^e^e  +  ^/lpP/lP  ‘^T,  +  (132) 

+  4>GcRcC  4>GC  +  2<l>,P„ftP  4>R/*  +  24>,Pe,T,  <^T,  +  2^ePe>r^  +  ^^ePetPM  ^PM 
+  2^ePetGC  <1>CC)  R 


/«  o\ 

where  R  =  (  )  (133) 

\0  R/ 


R 


r,r 


4>« 


(A  A  A  A  \  • 

,  rxv  .  rxv\.  r  .  r 

»*  X  r  )  r  =  V  =  ~ 

Ir  X  vl  Ir  X  vl  /  Irl  Irl 


(134) 


inertial  position  and  velocity  obtained  by  interpolating  off  of  the 
trajectory  at  tj 

6X6staU  transition  matrix^^'obtainedbyinterpoUUngoirof 
the  trajectory  at  tj 


partials  of  r  and  r  at  ty  with  respect  to  RP  interpolated  off 
of  the  trajectory 


partials  of  r  and  f  at  ^  with  respect  to>  T  computed  using 
equation  (79)  in  the  Range  Partial  Derivatives  subsec- 
tion  except  evaluated  at  ^  and  all  6  rows  are  required 


<1>PM  — 


partials  of  r  and  f  at  ^  with  respect  to  PM  computed  as  in 
equation  (80)  in  the  Range  Partial  Derivatives  subsec¬ 
tion  except  evaluated  at  tj  and  h  replaces  r  in  computing 
the  velocity  partials 
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4>cc  = 


partials  of  r  and  r  at  tj  with  respect  to  GC  interpolated  off 
of  the  trajectory 


Pet  PRPt  Pvt  PpMt  Pcct  PejiPt  Pe.Tt  PePMt  and  are  portions  of  the  full  covariance 
matrix  computed  by  equation  (125)  above.  Only  the  square  roots  of  the  diagonal 
elements  of  Prac  are  required  in  km  and  km/sec.  These  are  converted  to  meters  and 
mm/sec  for  plotting. 


The  corrections  and  standard  deviations  for  the  two  pole  coordinates  are  con¬ 
verted  to  km  at  the  Earth's  surface  by  multiplying  hy^anarth  and  then  to  meters  for 
plotting.  The  correction  and  standard  deviation  for  Kt  are  converted  to  msec/day. 


TOTAL  CLOCK  OFFSETS 

The  current  state  satellite  time  and  frequency  offsets  from  GPS  time  and  fre¬ 
quency  are  computed  by  adding  together  the  nominal  clock  offsets  and  the  current 
state  solved-for  clock  correction  at  tj,  i.e., 

~  +  (135) 

^irou,t(tj)  =  +  (135) 

The  nominal  clock  may  contain  jumps  and  is  computed  in  the  Filter  at  each  mini¬ 
batch  step.  The  total  satellite  clock  offsets  can  be  computed  in  both  (he  Filter  and 
Smoother.  The  SATRACK  project  requires  the  Smoother-derived  total  clock  offsets 
for  both  the  satellites  and  stations.  The  values  for  the  stations  are  obtained  using 
parallel  computations. 


RESIDUALS  AND  SIGNAL-TO  NOISE 


Residuals  after  fit  are  computed  in  the  Filter  for  each  mini-batch  interval  by 
linear  adjustment  of  the  original  residuals,  i.e., 


(O-C)adj.  =  (O-C)'  -A'^Xj 


(137) 


where  AX^  is  the  solution  for  all  states  and  A'  and  (O— O'  were  saved  before  being 
whitened  and  decorrelated.  The  si^al-to-noise  ratio  is  defined  as  the  square  root  of 
the  weighted  sum-of-squares  of  residuals  divided  by  the  number  of  observations.  For 
each  mini-batch  interval  it  is  computed  after  fit  by 


(S/N)j(0-C) 


Mj  J 


(138) 
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where  Oobs  is  the  actual  sigma  used  (may  be  different  from  the  input  value  because  of 
the  minimum  sigma  override)  and  Mj  is  the  total  number  of  observations  processed  in 
the  mini-batch  interval.  The  signal-to-noise  ratio  for  each  mini-batch  interval  in 
the  Filter  is  also  available  without  computing  adjusted  residuals  as  a  by-product  of 
the  square  root  information  implementation.  The  •  vector  in  the  information  array 
after  measurement  update  is  related  to  the  signal-to-noise  as  follows: 


(S/N  )j(Filter) 


(139) 


Residuals  and  signal-to-noise  after  fit  from  the  Smoother  for  each  mini-batch 
interval  are  obtained  by  reprocessing  the  same  observations  used  in  the  Filter.  This 
is  done  by  evaluating  the  observation  equations  using  both  the  nominal  and  solved- 
for  clock  information,  tropospheric  refraction  corrections,  updated  station  coordi¬ 
nates,  propagated  trajectories  (described  in  the  next  section),  and  updated  polar 
motion  information.  If  the  corresponding  Filter  execution  used  previously  computed 
total  clock  offsets  for  the  satellites,  this  same  information  is  used  in  generating 
residuals  from  the  Smoother.  Equation  ( 138)  is  used  with  (O  —  C)o4j.  replaced  by 
(O^C)smr.  to  get  the  Smoother  signal-to-noise  ratio  (S/N)j( Smoother)  for  each  mini¬ 
hatch  interval.  Then  the  overall  Smoother  signal-to-noise  ratio  is  computed  by 

r  N  n' 

S  Mj((S/N  )^(Smoother)] 


Overall  S/N  (Smoother)  = 


j  =  0 


(140) 


/  =  0 


where  N  -h  i  is  the  number  of  mini-batch  steps  (corresponds  to  to  thru  t/y).  The 
number,  mean,  standard  deviation,  and  RMS  of  residuals  for  the  entire  fit  span  by 
satellite,  station,  and  overall  for  each  data  type  are  computed.  Residuals  are 
converted  from  kilometers  to  meters  for  plotting.  Residuals  are  also  computed  for 
observations  not  processed  in  the  Filter  because  they  were  tagged,  did  not  pass  the 
elevation  am^le  tolerance  test,  or  were  deleted  by  pass  number.  On  option  residuals 
for  range  difference,  differenced  range,  and  doubly-differenced  range  data  can  also 
be  computed  if  range  data  for  the  same  stations  were  processed  in  the  Filter. 


TRAJECTORY  PROPAGATION 


All  orbit-related  parameter  corrections  are  transformed  into  inertial  position 
and  velocity  corrections  at  each  trigectory  timeline  usin^  linear  propagation  tech¬ 
niques.  These  corrections  are  added  to  the  reference  trajectory  positions  and  veloci¬ 
ties  at  each  timeline  to  produce  the  propagated  tr^ectory.  Earui-fixed  position  and 
veloc;  ty  are  obtained  by  transforming  these  improved  inertial  coordinates  using 
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improved  values  for  polar  motion  if  also  solved  for.  The  propagated  tre^ectory  span 
normally  corresponds  to  the  fit  span  with  at  least  4  extra  timelines  added  on  both 
ends  to  accomodate  the  interpolation  method.  However,  in  the  case  of  no  orbit- 
related  p  parameters  being  solved  for  (''batch'*  mode),  the  propagated  tregectory  is 
identical  in  len^  to  the  reference  tregectory  (except  it  has  fewer  extra  timelines  at 
each  end)  even  if  the  fit  span  is  a  subset  of  the  tr^ectory  span.  Also  in  this  case  the 
partial  derivatives  and  other  quantities  (primarily  solar  radiation  pressure  model 
related  items)  on  the  orinnal  reference  tregectory  can  be  copied  to  tine  propagated 
trajectory  if  required.  All  satellites  are  processed  simultaneously.  Improved  initial 
conditions  required  for  integrating  new  trigectories  can  be  computed  after 
propagation  is  completed  as  described  at  the  end  of  this  section. 

The  corrections  to  inertial  position  and  velocity  at  a  given  trsgectory  timeline 
Tt  are  computed  as  follows: 


aHTt) 

dHTt) 

dHTt) 

AKT,)  = 

- ak»  + 

AC,  -1- 

-  A«j 

(141) 
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where  the  orbit-related  parameter  corrections  at  to  are  used  if  Te  <  to,  the  corrections 
at  tj  are  used  if  i,  and  the  corrections  at  ts  are  used  if  Tf^tyv.  The  partial 

derivatives  of  position  and  velocity  with  respect  to  orbital  elements,  radiation  pres¬ 
sure,  thrusUs),  polar  motion,  and  gravity  coefGcients  are  the  same  as  those  defined 
in  the  SOLUTION  AND  DIAGNOSTICS  section  except  they  are  evaluated  at 
The  partial  derivatives  of  position  and  velocity  with  respect  to  stochastic  radiation 
pressure  and  ^avitational  acceleration  parameters  are  the  same  as  those  defined  in 
the  STATE  EQUATIONS  section  with  j  replaced  by  Tt  ( equations  (22)  and  (23) 
for  Kr  and  equations  (29)  and  (30)  for  C).  These  partials  are  zero  if  7^ =t>  This  is  a 
result  of  the  pseudoepoch  state  formulation  of  the  equations.  If  a  Filter  ^atch"  mode 
propagated  trajectory  is  being  created  the  orbital  element  and  orbit-related  y  param¬ 
eter  corrections  at  ts  are  usea  at  all  tngectory  timelines.  If  a  Smoother  propa^ted 
trsgectory  is  being  created  the  orbit-related  y  parameter  corrections  at  ts  are  also 
usM  since  no  smoothing  of  these  corrections  is  possible. 

The  improved  inertial  position  and  velocity  at  a  given  tngectory  timeline  are 
then  given  by 


n^pro^^Tt)  =  r/:,/^^(T,)-fAr<T,)  (143) 

hmprootiATt)  —  ^R»f.(Tt)  +  ^Tt)  (144) 
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The  improved  Earth-fixed  position  and  velocity  are  then  given  by 

ricfiTf)  =  ABCD(Tt)r,^pro.aTt)  il46) 

hFiTt)  =  ABCI>{Tt)r,n,provt<^Tt)-k-AWBCD(Tt)r,^^root^Tt)  (146) 


where 


/  0  w  0\ 

W  =  -w  0  0  ) 

\  0  0  0/ 


(147) 


When  propagating  a  trajectory  based  on  the  Filter  corrections  the  improved  polar 
motion  information  is  used.  For  p  and  q  the  corrections  at  tj  are  used  at  trajectory 
timelines  such  that  +  For  A*(  the  same  interval  applies  and  the  At 

correction  is  converted  to  a  correction  to  At  by  multiplying  it  by  Tf-  to. 


Improved  initial  conditions  are  computed  if  required  after  trajectory  propaga¬ 
tion  is  completed.  For  each  mode  of  operation  the  computations  are  slightly  differ¬ 
ent.  In  all  modes  the  improved  inertial  position  and  velocity  are  obtained  by  inter¬ 
polating  off  of  the  propagated  trajecto^  at  the  selected  times.  In  the  "batch^'  mode 
the  ARP\  values  are  added  to  the  nominal  radiation  pressure  parameter  values  to 
get  the  improved  radiation  pressure  parameter  values-the  same  for  all  initial 
condition  times.  A  similar  computation  is  done  for  thrust  parameters.  In  the 
"Filter"  mode  both  the  ana  ARP.  values  just  before  the  time  for  initial  condi¬ 
tions  are  added  to  the  nominal  radiation  pressure  parameter  values  to  get  improved 
values.  A  similar  computation  is  done  for  thrust  parameters.  In  the  "Smoother" 


N 


mode  the  ARP/^  values  and  the  average  AICm,  values 


/=0 


are  both  added 


N  +  J 


to  the  nominal  radiation  pressure  parameter  values  to  get  improved  values.  These 
radiation  pressure  parameter  values  would  normally  be  used  to  predict  a  reference 
trsjectory  for  a  future  span.  This  is  why  the  average  stochastic  radiation  pressure 
corrections  are  used.  The  thrust  corrections  are  a^ed  to  the  nominal  thrust  values 
to  get  improved  values.  Also  the  gravity  coefficient  parameter  corrections  are  added 
to  their  nominal  values  to  get  improved  values. 


SATRACK  COVARIANCE  MATRIX  GENERATION 


Special  covariance  matrices  are  required  for  the  SATRACK  application  of  the 
MSF/S  system  of  programs.  These  matrices  relate  8  parameters  (position,  velocity, 
time  offset,  and  frequency  offset)  for  each  satellite  to  the  same  parameters  for  every 
satellite  at  a  given  time  ( intersatellite  covariances)  and  at  up  to  three  different 
times  (intertime  covariances).  The  full  covariance  matrix  required  is  structured  as 
follows: 
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PSATHACK 


8Nsy  SNjfy  BNgy 


(148) 


where  the  subscripts  refer  to  times  and  Tj  and  Nsy  is  the  number  of  satellites. 

Pn,  P22y  and  P33  are  symmetric  submatices  and  P2;  =  Pfa,  P^i  =  Pfj,  and  Pj2=  P23. 
Each  submatrix  is  further  divided  into  8X8  blocks  where  blocks  along  the  diagonal 
relate  a  satellite  to  itself  either  at  the  same  time  or  at  two  different  times.  Eacn  8X8 
block  above  or  below  the  diagonal  blocks  in  each  submatrix  relate  one  satellite  to 
another  satellite  either  at  the  same  time  or  at  two  different  times.  The  diagonal 
submatrices  Pu,  P22,  and  P33  are  computed  exactly  but  the  off-diagonal  submatrices 
Pt2,  Pt3y  and  P23  are  approximated  in  such  a  way  that  they  would  be  exact  if  the 
process  noise  terms  were  zero. 

As  indicated  in  Figure  3  in  the  TIMELINE  DEFINITIONS  section,  the  times 
of  interest  T/,  and  T3  may  not  be  exactly  on  mini-batch  steps.  Each  is 

/  At'***' 

associated  with  the  mini-batch  interval  such  that  T^e  I  tj  —  --  .  tj  +  J  . 

Let  R  m  =  1,2,3  denote  the  Pj  'matrices  for  these  mini-batch  times.  Each 
matrix  is  of  dimension  Nt’otXNt’o?'  and  is  upper  triangular.  For  each  time  7,„  a 
state  transition  matrix  ^(T^)  is  defined  as  given  below.  This  matrix  maps  the 
parameter  corrections  at  into  corrections  in  position,  velocity,  satellite  time  offset, 
and  satellite  frequency  onset  at  Tn. 


G 


0 


corresponds 

toCfi 


0 

♦ 

1 

corresponds 

Cms 


0 


RP 


T  iPjlf 


^)i 


GCnSNsy  (149) 


corresponds 

toS 


Columns  are  present  only  for  the  parameters  included  in  A  '(7^). 
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(150) 


G  and  RP  have  the  same  size  and  structure  as  Kh 


K/i^  =  partials  of  position  and  velocity  at  with  respect  to  stochastic 

radiation  messure  parameters  at  t,  given  by  equations  (22)  and  (23)  in 
the  STATE  EQUATIONS  section  with  tj+ 1  replaced  by 

Gi  =  partials  of  position  and  velocity  at  with  respect  to  gravitational 
acceleration  parameters  at  tj  given  by  equations  (29)  and  (30)  in  the 
STATE  EQUATIONS  section  with  tj+t  replaced  by 

RPi  -  partials  of  position  and  velocity  at  tj  with  respect  to  radiation  pressure 
parameters  obtained  by  interpolating  off  of  the  trigectory  at  T„ 
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=  partials  of  position  and  velocity  at  T„  with  respect  to  orbital  elements 
obtained  by  interpolating  off  of  the  trigectory  at 


The  rows  for  each  Tj  are  determined  by  which  satellites  have  thrusts. 

Ti  =  partials  of  position  and  velocity  at  wi  respect  to  thrust  computed 
using  equation  (79)  in  the  Ryge  Partial  Derivatives  subsection 
except  evaluated  at  and  all  6  rows  are  required 


59 


NSWCTR  87-187 


NpM 


0 

PM  =  PMj 

SNsvXNpM  0 


GC  has  the  same  size  and  structure  as  PM 


(156) 


PMi  =  partials  of  position  and  velocity  at  with  respect  to  polar  motion 

computed  as  in  equation  (80)  in  the  Range  Partial  Periyatiyes  sub¬ 
section  except  evaluated  at  T„  and  f  replaces  r  in  computing  ue  velocity 
partials 

GCi  =  partials  of  position  and  velocity  at  with  respect  («  gravity 
coefficients  obtained  by  interpolating  off  of  the  treoectory  at  T„ 


The  following  product  matrices  are  then  formed: 

S(TJ  =  <t>(T„)R'(Tj  m  =  1,2,3  (156) 

SSsv^f^TOT  SNsv^Ntot  NroT^N'ror 

Then  the  covariance  submatrices  of  Psatrack  are  given  as  follows: 

Pit  =  S(T,)S{Tjr  (157) 

Pi2  =  S(T,)S(T2)^  (158) 

Pt3  =  S(T,)S(T3f  (159) 

P22  =  S(T2)S(T2)^  (160) 

P23  =  S(T2)S(T3y  (161) 

P33  =  S(T3)S(T3V  (162) 


Each  is  an  8Nsv  XSNsv  matrix  where  the  eight  parameters  are  position,  velocity, 
satellite  time  offset,  and  satellite  frequency  offset  in  units  of  km,  km/sec,  psec,  and 
ppm  respectively.  These  matrices  are  then  scaled  to  be  in  units  of  m,  m/sec,  sec,  and 
sec/sec. 
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APPENDIX  A 

STATE  EQUATION  SIMPUPICATIONS 
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The  puriMse  of  this  ap^ndiz  is  to  describe  the  assumptions  and  definitions 
made  in  adopting  the  specisilized  form  of  the  state  equations  for  the  orbit-related 
parameters  in  the  MSF/S  system.  The  Hrst  simplification  involves  assumptions 
about  the  white  noise  sources  driving  the  estimated  orbit-related  states  and  affectd 
the  form  of  the  process  noise  covariance  matrix.  The  second  simplification  involves 
the  deHnition  of  the  pseudoepoch  state  variables  and  affects  boUi  the  state  transition 
matrices  required  for  the  Filter  propa^tion  step  and  the  observational  equation 
partial  derivatives  required  for  theFilter  measurement  update  step. 

To  simplify  this  discussion,  consider  the  set  of  state  equations  for  a  single  satel¬ 
lite  that  includes  current  state  TOsition  and  velocity  parameters,  x',  and  one  stochas¬ 
tic  orbit-related  parameter,  p.  The  system  of  stochastic  differential  equations 
describing  corrections  to  these  parameters  is  given  by 


where  b  = - ,i.e.,Sp  is  modeled  as  a  first-order  Gauss-Markov 

^  process  (see  Appendix  B) 


ax' 

Fp  =  —  =  partial  derivatives  of  velocity  and  acceleration  with 
dp  respect  to  p 


ax' 

f;  =  — =  partialderivativesofvelocity  and  acceleration  with 
rix'  respect  to  position  and  velocity 

and  it  is  assumed  th^  t  changes  in  x'  cause  negligible  changes  in  p.  Also  each  com- 

ponent  is  a  white  noise  process.  Assume  that  =0,  i.e.,  there  is  no  white  noise  driv¬ 
ing  the  8x'  states  directly.  Then  the  discrete  system  equivalent  to  equation  (Al)  is 
given  by 


(tj+t-tj) 


where  M  =  e  ^  (see  Appendix  B) 

5x'(tj>i) 

~  — - =  partial  derivatives  of  position  and 

dpiXj)  velocity  at  tj+i  with  respect  to  p  at  tj 

dx‘{tj+j) 

y'xi  =  Kitj+htj)  =  -  =  partial  derivatives  of  position  and 

ax'(t^)  velocitv  at  tj  /  with  respect  to  position 
and  velocity  at  tj 
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(A-tj) 


v;(A,tj)  v;(A,tj) 


a-tj) 


rp(X,tj)  w, 


Therefore,  W2  is  a  non-zero  vector  as  a  result  of  the  discretization  process  even  thoucrh 

A  «srW « «%  M  A  m  m  ^  ^  A.  ^  ^  i  x1__ A.* _ _ .j  •  « 


smoothable  since  they  are  dynamically  related  to  the  Ap  state  through  Ae  Vp 
matrix.  This  assumption  results  in  the  absence  of  both  a  process  noise  covariance 
matrix  for  the  six  orbit  states  and  a  process  noise  cross-covariance  matrix  between 
these  six  states  and  the  stochastic  orbit-related  states.  Since  the  orbit  states  are  non¬ 
stochastic  (not  p  parameters),  a  considerable  savings  in  array  storage  results.  This  is 
the  case  because  the  information  array  required  for  the  Filter  propagation  step 
(equation  (111))  includes  two  rows  ana  columns  for  each  p  parameter.  This  reduction 
therefore  allows  more  satellites  to  be  processed  simultaneously. 

Expanding  the  state  equations  given  in  equation  (A2)  to  include  all  orbit- 
related  stochastic  states,  Ap,  and  orbit-related  bias  states,  Ay,  and  again  assuming 
^2  =  0,  gives  the  following  state  equations: 


Ax'  = 


V'  v;  v: 

Pj  *j  Jj 


where  M 


=  diag\e  X|  / 


and  Vv  = 


yp(tj+i,tj)  =  partial  derivatives  of  position  and  velocity  at  t^+j 
with  respect  to  p  at  tj 

=  partial  derivatives  ofposition  and  velocity  at  tj+/ 
with  respect  to  position  and  velocity  at  tj 

Vy{tj+i,tj)  =  partial  derivatives  of  position  and  velocity  at  tj+i 
with  respect  to  y  at  tj. 
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The  subset  of  these  equations  involving  the  Ax' states  is  given  by 


=v;  Ap(tj)  +  v;(tj+,,tj)Ax'(tj)  +  v;(tj>i,t^)  Ay(ty)  (A5) 

Define  the  pseudoepoch  state  variables  Ax^  by 

Ax-ftj)  =  V,  (tj.  To)  Ax^ + V,  {tj.  To)  Ay(tj)  ( A6) 

where  V,  (tj,To)  =  partial  derivatives  of  position  and  velocity  '\t  tj  with  respect  to 
either  orbital  elements  or  position  and  velocity  at  Uie  trajectory 
epoch  To 

and  Vy  {tj,T o)  —  partial  derivatives  of  position  and  velocity  at  tj  with  respect  to 

the  orbit-related  bias  parameters  at  To. 

Substituting  for  Ax'fty)  and  Ax‘(tj+i)  from  equation  (A6)  into  equation  (A5)  results  in 
the  following  equation: 

V*(tj+,,To)  Axj+i  +  V,  (tj+i,To)  Ay(t^+,)=  (A7) 

Ap(tj)-I- Vi  (tj+j,  t;)V,(t^,To)Axj-l- tj)Vy(tjJo)Ay(tj)  +  V^  (tj+j,  tj)Ay(tj) 

Multiplying  both  sides  of  equation  (A7)  by  Vxitj+j,  To)  and  rearranging  terms 
results  in  the  following  equation: 

Ax;+,  =  V,'(tj+,,To)v;(t;+,,t,)V,(t^,To)Ax;  +  (tj+ j>  To)Vp  (tj+ j,tj)Ap(tj)  (A8) 

+  V;'(tj+;,To)[(v;(tj+,,t;)V,(tj.To)-^v;(t,+,,t,))Ay(tj)-y,(tj+,Jo)Ay(t,^^ 

Based  on  the  properties  of  state  transition  matrices,  the  following  identity  exists: 

v;(t;+/,t,)V,(tj,To)  =  V,  itj^iJo)  (A9) 

Therefore  the  matrix  multiplying  Ax,  in  actuation  (A8)  reduces  to  an  identity  matrix. 
Also,  based  on  the  properties  of  state  transition  matrices  and  Ay(tj+i)-  Ay  (t^)  from 
equation  (A4),  the  lollowing  identity  exists: 


[Vx(tj+I ,tj)Vy{tj,To)  +  Vyitj^t ,tj)]  Ay(V  =  Vy  {tj+i,To)Ay{tj+i) 


(AlO) 


Therefore  the  last  term  in  equation  (A8)  vanishes  and  equation  (A8)  reduces  to: 

Axj+i  =  Axj+Vxitj^,,To)ypitj+t,tj)Ap{tj)  (All) 

=  Axj+VpAp(tj) 

where  Vp^  =  Vx(tj+j,To)Vp(tj+i,tj)  (A12) 
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Therefore  the  state  equations  given  by  equation  (A4)  and  the  pseudoepoch  state 
variables  definition  ^ven  by  equation  (A6)  result  in  the  following  simplified  state 
equations  for  orbit-related  parameters: 

/Ap\  /M  0  0\  /Ap\  /wA 

Ax  =  /  0  j  Ax  j  -H  j  0  j  (A13) 

\Ay  /  \0  0  /  /  \Ay  /  \o  / 

To  be  consistent  with  these  equations,  the  observational  equation  partial 
derivatives  required  in  the  Filter  measurement  update  step  must  be  computed  with 
respect  to  these  variables.  For  a  given  observation  at  time  tob«,  the  partials  with 
respect  to  the  orbit-related  p  parameters  involve  the  Vj,  Uob„  tj)  matrix,  the  partials 
witn  respect  to  the  pseudoepoch  state  orbit  parameters  involve  the  V,  (tob».  To) 
matrix,  and  the  partials  with  respect  to  the  orbit-related  bias  parameters  involve  the 
Vy  (tob„  To)  matrix. 


A-6 


NSWCTR  87-187 


APPENDIX  B 

FIRST-ORDER  GAUSS-MARKOV  PROCESSES 
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The  stochastic  differential  equation  deftning  a  continuous  first-order  Gauss- 
Markov  process  is  given  by 

8p  =  — —  8p+w  (Bl) 

T 

where  t  =  f’ecorrelation  time  (seconds) 

w  =  white  Gaussian  noise  of  mean  zero  and  spectral  density  q,  i.e.,  £(w) = 0 
and  E[w(t)w(s)]=q8(t— s). 

The  continuous  linear  variance  differential  equation  corresponding  to  this  stochastic 
process  is  given  by 

a*  =  — —  o*-hq  (B2) 

X 

Setting  the  variance  rate,  6^,  to  zero  and  solving  for  a*  gives  the  steady-state  variance 
for  this  random  process  as 

a*  =  £(8p(t)*)  -\q  (B3) 

m 


Also  the  mean  value,  p,  of  this  process  is  zero,  i.e., 

p  =  e(8p(t))  =  0  (B4) 

The  one-sided  power  spectral  density  for  a  Gauss-Markov  process  is  given  by 

4,(u)  =  (B5) 

<a*  +  P» 


where  P  =  —  and  u  =  2nf. 

X 

Its  corresponding  autocorrelation  function  is  given  by 

<l>(t)  =  a*e-P'*'  (B6) 

A  Gauss-Markov  process  is  used  to  approximate  a  band-limited  process  with  a  flat 
spectral  density  over  Uiis  bandwidth. 

The  discrete  equivalent  of  this  continuous  process  is  gpven  by  the  first-order 
stochastic  difference  equation 

Apj+i  =  4>Ap,  +  wy  (B7) 

where  4>  =  state  transition  matrix 

{w,}  =  white  Gaussian  noise  sequence  of  mean  zero  and  variance  Q,  i.e., 

£(wy) = 0  and  EiyijMfk) = Q8(/-*). 
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<|>  satisnes  the  following  differential  equation  and  initial  condition: 

♦  =  -i* 

t 

4»(W  =  1 

The  solution  to  this  differential  equation  is  simply 

4>  =  e 


X  =  •  t 

The  process  noise  variance,  Q,  is  then  given  by 


Q  = 


tj+i 


rt>+i  2(X—tj) 

s  t 


dA 


(B8) 

(B9) 


(BIO) 


2(k-tj) 


**  2* 


(BID 


2At 

0*  [l-e  X 

For  x=<B,  the  stochastic  differential  equation  (Bl)  reduces  to  that  for  a  random 
walk  process,  i.e.. 


6p  =  w 

The  discrete  equivalent  is  then  g^ven  by 

Apj+j  =  Apj-hvj 

For  this  case  Q  is  derived  starting  with  equation  (Bll)  as  follows: 


(B12) 


(B13) 


Q  =  Urn  q  — 
2 


_2At' 
1-e  X 


(B14) 
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=  q  Urn  (  £  conjtanti/—  ]  ) 

1  =  1  \x/ J 


=  qAt 

If,  in  addition,  q=0,  then  the  stochastic  differential  equation  (B12)  reduces  to  that 
for  a  random  constant,  i.e.. 


8p  =  0 

The  discrete  equivalent  is  then  given  by 
^Pj+i  =  Apj 

and  Q  =  0 


(B15) 


(B16) 

(B17) 
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APPENDIX  C 

APPROXIMATION  OP  PARTIAL 
DERIVATIVES  OF  POSITION  AND  VELOCITY 
WITH  RESPECT  TO  STOCHASTIC  ORBIT-REIJVTED  STATES 
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Hie  partial  derivativeai^poaitioa  and  vebeity  at  t  with  reepeet  to  the  stochaa- 
tic  radiation  preiaure  and  graeitaticMial  acceleration  parameters  at  ti  are  approxi¬ 
mated  using  a  Taylor  series  expansion  method.  These  partials  only  nave  to  be 
accurate  for  values  of  t  s  f ,  i.e.,  for  a  time  ^an  of  at  most  one  mini-batch  interval 

(usually  s  1  hour).  Using  p  to  di^gnata  the  stochastic  parameters  (either  ig  or  G) 
and  r  ^d  f  to  designate  pomtion  and  velocity  reepectivmy,  the  aysUm  of  dififorential 
equations  defining  the  corrections  to  the  parameters  (with  the  stochastic  terms  set  to 
sero)  used  to  derive  the  state  transition  matrix  containing  the  required  partial 
derivatives  is  given  by 


ar 

Ft  =  —  —  partial  derivatives  of  acceleration  with  respect  to  p 

dp 


and  C  =  partial  derivatives  ofacceleration  with  respect  tor 


It  is  assumed  that  changes  in  r  and  f  cause  negligible  changes  in  p  and  that 
acceleration  is  not  a  function  of  velocity. 

Define 


F 


where  Fg 


and  Fg 


(C2) 

(C3) 

(C4) 


The  state  transition  matrix,  is  the  solution  of  the  following  system  of 
differential  equations  and  initial  conditions: 


4>(t,ty)  =  F(t)<^t,ty) 


(C5) 


=  / 


(C6) 
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Approximating  ^t,ty)  by  a  second-order  Taylor  series  about  tj  and  assuming  F = 0, 
i.e.,  F  is  constant  for  short  time  intervals,  results  in  the  following: 


(t-t/ 


(C7) 


r  1 

=  1+  F(t^)<|Ktj,tjKt  -  tj)  +  [F(t^)<HM>) + J  — — 


(t-t/ 


/ 


B* 


0 


where  F*(tj)  = 

\F3(tj)Faitj)  +  FAtj)B  Fiitj)^ 

Partitioning  ({k  the  same  way  that  F  was  partitioned  in  equation  (C2)  gives 

'♦i(t.tj)  0 


<|Kt.tj)  =1 


Then  it—tjY 

A  l+fl(t-tj)  +  B* - 


(C8) 


(C9) 


(CIO) 


However,  in  the  limit  as  the  number  of  terms  increases,  the  right-hand  side  of  this 

expression  converges  to  diag\e  x,  Iso  no  approximation  is  necessary.  Also, 

-  (t-t.)* 

(Cll) 


*  /+Fa(t;Mt-tj)+ Fa(tj)  - 


This  anproximation  is  also  not  required  since  ^j(t,  t^)  is  just  the  partials  of  portion 
and  velocity  at  t  with  respect  to  position  and  velocity  at  ty^whicn  can  be  obtained 
exactly  by  proper  manipulation  of  partials  interpolated  on  of  the  trajectory. 
Therefore,  the  only  submatrix  of  ^(t,ty)  requiring  approximation  is  and  it  i 
given  by 


IS 


4>2(t,t^) 


(' 


F3(tj)Faitj)  +  F^tj)Bj 


(t-t/ 

~2 


(C12) 
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Therefore 


Mtj)  it-tjf 

-  *  Ff(ty)  - 

ap(v)  * 

dHtj) 


(t-tjY 


(C13) 


(C14) 


aRt,) 


For  the  radiation  pressure  (Ka)  parameters,  F|(tp=  - is  computed  as 

follows:  dKi^tj) 


ritj)  =  R, 


iO  '^ shape  cos 
Kr^  10  '‘Shape  sin  Kr^ 


(C15) 


aRt,) 

aKA(t>) 


/a,  lO’^shaptcosKR^  ^KR^10  ’*3kap$  sin  Kr^ 
=  R,  0  10 shape  sin  Ka,  Kr^10  ‘*  shops  cos  Kr, 


(C16) 


where  R,  =  matrix  required  to  transform  between  the  body-axis  and  inertial 
Cartesian  reference  systems  obtained  from  the  tn^ectory  at  tj 


a,-KA,10  “shape  cos  Ka, 


M 

.  O,  = 


(C17) 


=  inertial  accelerations  due  to  the  radiation  pressure  model  only  and  not 
including  y-axis  and  Ka,  contributions,  given  in  the  body-axis  x  and  z 
directions  respectively  at  tj, 

[  a^  ]  =  inertial  acceleration  at  tj  due  to  radiation  pressure  in  the  body-axis 
\a«/  directions  obtained  from  the  trajectory 

Kr  =  nominal  radiation  pressure  parameter  values  from  the  trajectory 

shops  =  fraction  of  the  sun’s  disk  unobstructed  by  any  eclipsing  body  (Earth, 
Moon,  or  both)  obtained  from  Uie  traject^  at  tj 

aRt^) 

For  the  gravitational  acceleration  (G)  parameters,  F/itj)  = - is  computed 

as  follows:  dG(tj) 


Rtj)  —  Rrac^ 


(C18) 
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ar(t,) 

-  =  8jiac 


where  Rhac 


/ 


f  X  V 

r — 
Ir  X  vl 


rx  V  \ 

Ifx^/ 


(C19) 

(C20) 


=:  matrix  required  to  transform  between  the  RAC  and  inertial 

Cartesian  reference  frames  at  tj 

r,v  =  position  and  velocity  at  interpolated  off  of  the  trs^tory 
denotes  a  unit  vector) 
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APPENDIX  D 
CLOCK  MODELS 
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The  current  state  satellite  clock  offsets  are  modeled  as  a  third-order  linear 
system  with  white  noise  inputs  given  by  the  following  stochastic  differential 
equations: 


where  6V 


0\  /Sir* 


frequency  drift 
frequency  offset 
time  offset 

white  noise  process  of  mean  zero  and  spectral  density  qi,  i.e., 
£(  wy  (t))  =  0  and  Etw?  (t)wy  (s))  =  q.  8(t  -  s),  i  =  1,2,3 


The  discrete  equivalent  of  tliis  continuous  model  is  given  by 

/A*;\  /I  0  0\  /A*t*\  /wi\ 

'  a;  I  =  (  At  J  0  I  I  A?  I  +  (  wi  I 

V-i,  -  P  w  vw 


(D2) 


where  At 

{w;} 


tj+i-'tj 

white  noise  vector  sequence  of  mean  zero  and  covariance 
matrix  Q  derived  as  follows: 


/qi  0  0 

I  0  02  0 


1  X  ^ 
2 


dk  (D4) 


Ql 

qjA 

q,A 

q,A*  +02 

A* 

■"i 

A*’ 

qj  -  +02X 
2 

Qi — dA 

/  X..  / 

Oi  --  -rq^A^+qj J 
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q,At 


1-7 


it’  it’  , 

2  3 


At'  At^  .  At^' 

Ql —  Ql—  +  Q2  — 
6  8  2 


In  terms  of  the  pseudoepoch  state  clock  offsets  actually  implemented  in 
MSF/S  system,  the  model  is  ^ven  hy 


where  the  subscript  o  denotes  a  correspondence  with  the  Ht  span  epoch,  to,  an 
{wj}  =  white  noise  vector  sequence  of  mean  zero  and  covariance  matrix  Qo 
from  Q  as  follows: 


Then  Arj+i  =  4>j+/  Ato^^,  =  <J)y+/ Ato^  +  ^j+i  Wj 

The  state  transition  matrix  in  equation  (D2)  is  equivalent  to  ^j+}^Jso  that 
equation  (D2)  can  be  rewritten  as 

ATj+i  =  +  +  W} 

Combining  equations  (D6)  and  (D7)  results  in  the  following  equality: 

4>j  +  / Aro^+  ^j+t  W/  =  +  wj 

Multiplying  each  side  by  its  transpose  gives 

4>J+/  +  ^J  +  1  4‘J+J  +  ^j  +  t  wJ  -I-  ^j+t  4>J+J  = 

4>J  +  i + »/  wj  ^ wj^+  wj  ^tJ  +  / 


D-4 


the 


(D5) 


derived 


(D6) 

(D7) 

(D8) 

(D9) 
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Taking  expected  values  and  replacing 
AToj)  by  Po^, 

E(wjwj)  by  Qo^, 

E{Atow])  ,  Einfj^xo),  E(Atj  Wj*'),  and  £(wj  ArJ)byO, 

E(^rj  AtJ)  by  Pj , 
and  E(wjwJ^  by  Q  gives 

<Pj+ /  ^Oj  4>/+i + ^j+iQoj  4>J+  /  =  4»j + 1  ^jPj  4>/  1  +  Q  (DIO) 

Multiplying  each  side  of  this  equation  on  the  left  by  i  and  on  the  right  by  4>j>i  and 

substituting  for  gives 

Poj+Qo^  =  Po+  (Dll) 


Therefore  Qoj  =  ^nj'+iQ^V+i  (D12) 

The  square  root  of  Qo^  is  actually  required  hy  the  Filter  algorithm.  The  upper 
triangular  square  root,  P,  of  Q  '  such  that 

Q  '  =  R^R  or  Q  =  R  '  R  ''  (D13) 

is  first  com^ted  using  Cholesky  decomposition  (see  Appendix  E).  Then  combining 
equations  (D12)  and  (D13)  gives 

Qo,=  $,>/«' (D14) 
Inverting  gives 

(D15) 


Therefore  the  square  root  of  Qo  is  just  R  x  where  R  is  invariant  (a  function  of  At 
only)  and  ^j+j  is  lower  trian^lar.  The  product  R^j+ ;  is  a  full  matrix  and  must  be 
upper  triangularized  usinjr  Householder  transformations  (see  Appendix  G)  before 
being  used  in  the  propagation  array. 

To  be  consistent  with  the  clock  model  states  as  given  in  equation  (D5),  the 
observational  partial  derivatives  required  in  the  Filter  measurement  uj^ate  step 
must  be  computed  with  respect  to  these  pseudoepoch  state  parameters.  For  a  given 
observation  at  time  to6$,  the  partials  must  thererore  involve  the  third  row  of  the 
matrix  with  tj  replaced  by  tob»- 
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The  current  state  station  clock  offsets  are  modeled  as  a  second-order  linear 
system  with  white  noise  inputs,  i.e.,  the  same  as  the  satellite  clock  model  except  that 
all  SV-related  terms  are  absent.  For  the  discrete  version  of  this  model  the  Q  matrix  is 
given  by 


(D16) 


where  qt  =  white  noise  spectral  density  for  the  continuous  frequency  offset  state,  5r 
and  02  =  whi  te  noise  spectral  density  for  the  continuous  time  offset  state,  8r. 

The  matrix  is  given  by 

fc  =  (  *  (D17) 

\tj-to  IJ 

With  these  definitions,  all  of  the  discussion  of  the  satellite  clock  model  above  applies 
to  the  station  clock  model  also. 


Assuming  thatQi  =0  for  the  satellite  clock  model,  i.e.,  frequency  drift  is 
modeled  as  a  random  constant  over  the  entire  fit  span,  the  model  for  me  frequency 
offset  state,  Ai,  is  equivalent  to  a  constant  plus  the  integral  of  the  frequency  drift 
state,  Af,  plus  the  integral  of  white  noise.  The  model  for  the  time  offset  state,  Ar,  is 
then  equivalent  to  a  constant  plus  the  integral  of  the  frequency  offset  state,  Af,  plus 
the  integral  of  white  noise.  Therefore,  the  noise  terms  integrated  in  determining  the 
time  offset  state  are  integrated  white  noise  (random  walk  frequency  noise)  and  white 
frequency  noise. 


The  spectral  densities  of  these  noise  terms  can  be  directly  related  to  the  Allan 
variance  wnich  is  used  to  characterize  the  statistical  frequency  fluctuations  of 

atomic  and  crystal  clocks.  The  Allan  variance,o*(r),  is  defined  by 

o^Ct)  =  -  E[i^k*iU)-9kU)Y]  (D18) 

2 

fU  +  t 

where  ~  \  y(t)dt  =  average  fractional  frequency  over  i  seconds  (D19) 

■  i. 


and  yit)  =  -  =  instantaneous  fractional  frequency  (D20) 

2nVo 

=  rate  of  change  of  phase  divided  by  2n  times  the 

nominal  frequency  v^. 


D-6 


NSWCTR  87-187 


The  Allan  variance  corresponding  to  the  model  with  qi  =  0  is  given  hy  the  sum  of  the 
individual  Allan  variances  for  each  noise  source  as  follows: 


(D21) 


The  Hrst  term  corresponds  to  the  white  frequency  noise  and  the  second  term 
corresi^nds  to  the  random  walk  frequency  noise.  This  is  depicted  in  Figure  Dl, 
which  is  a  typical  plot  of  the  square  root  of  the  Allan  variance.  Flicker  noise,  which 
corresponds  to  a  horizontal  line  on  this  plot,  cannot  be  represented  exactly  by  this 
model.  However,  the  white  noise  term  (with  spectral  density  q^)  can  be  chosen  to 
exactly  match  the  left-hand  portion  of  a  theoretical  Allan  variance  curve.  Then  the 
q2  spectral  density  can  be  selected  optimally  so  that  the  minimum  point  of  the 
combined  curve  lies  on  the  flicker  noise  portion  of  the  theoretical  curve. 


_ FLICKER 

NOISE 


FIGURE  Dl.  CORRESPONDENCE  BETWEEN  CLOCK  MODEL  SPECTRAL  DENSITIES 

AND  ALLAN  VARIANCE 
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APPENDIX  E 

CHOLESKY  DECOMPOSITION 


E-1 
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Lower  triangular  Cholesky  decomposition  is  used  to  compute  the  square  root  of 
the  inverse  ofthe  process  noise  covariance  matrix,  Q,  for  each  clock  modd,i.e.,  Q  ’  is 
factored  into  the  form  R^R  where  R  is  uj^per  trian{^lar.  The  algorithm  actually  com- 

?utes  a  lower  triangular  matrix  L = R^  i.e.,  Q  ‘  =  LL**.  Let  q-fiQ  'in  X  n)  and 
ij  eL(n  Xn),  n = 2  or  3.  Then  L  is  computed  as  follows  starting  with  /=  1: 

For  colunmj  define 

fjj  =  (q,v)'  (El) 

If  7 =n,  the  procedure  is  complete. 

Otherwise,  for  each  subsequent  row  k  define 

^kj  =  q*y^jV  ^  ~  /+l...,n  (E2) 

This  completes  the  definition  of  this  column  of  L  and  the  Q  '  matrix  is  then 
adjusted  as  follows: 

For  all  subsequent  columns  and  all  subsequent  rows  define 

Qijt  ~  k  —  7'^ l...,n  i  —  (£3) 

Then  go  to  the  next  colunm,  /+ 1,  and  repeat  this  procedure. 
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APPENDIX  F 

WHITENING  AND  DECORRELATION  OF  RANGE 
DIFFERENCE  OBSERVATIONS 


F-1 
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Range  difference  observations  are  derived  bv  one  of  two  measurement  tech¬ 
niques,  both  of  which  result  in  pairwise  correlated  observations.  The  first  technique 
involves  differencing  two  consecutive  accumulated  Doppler  counts  provided  that  the 
count  was  continuous,  i.e.,  no  losses  of  lock  or  cycle  slips  occurred  during  this  inter¬ 
val.  The  second  technique  involves  differencing  two  consecutive  pseudorange  meas¬ 
urements.  In  either  case,  two  consecutive  ranm  ^fference  observations  (assuming 
no  losses  of  lock  or  cycle  dips  for  the  Doppler-derived  type)  have  one  measurement  in 
common  and  are  therefore  correlated. 

Let  RD  represent  a  measurement-noise-free  range  difference  observation  and 
let  ADR  represent  a  noise-free  Doppler  count  type  measurement  or  pseudorange 
measurement.  Then  from  three  consecutive  ADR  measurements,  two  RD 
observations  can  be  derived  as  follows: 

RDi  -hvi  =  {ADRt  +vj )  -  (ADRo+vH )  (FI) 

RDs-l-vi  =  (ADRa-l-v^)  -  (ADR,-l-v?)  (P2) 

where  v"  represents  the  zero  mean  white  measurement  noise  of  varianceoAn/t,  on  the 
ADR  measurements  and  vj  represents  Uie  resulting  zero  mean  measurement  noise  on 
the  RD  observations.  The  RO  observation  measurement  noise  variances  are  derived 
by  squaring  each  side  of  equations  (FI)  and  (F2)  and  taking  expected  values  to  get 
the  allowing: 

Ofto,  =  E(vi*)=  E(vj^+  B(vo*)=  OadRj+  Oadr^  (F3) 


orOj  =  E(va^=  £(v2^+  E(vj‘)=  0/|/>r,+  Oadr,  (F4) 

The  measurement  noise  covariance  is  derived  by  multiplying  equation  (F2)  by 
equation  (FI)  and  taking  expected  values  to  get  the  following: 

^RDiJiDg  —  tXvivj)  =  —  E(v/^  =  —Oadr,  (F5) 

Then  the  correlation  coefficient  is  given  by 


pRDi.RD,  - 


_  OgD,.  RD^  _ 


i 

Oadr, 


(Ord,  Ord^)*  Ord,  Ofto, 


(F6) 


Let  Oi = ord^  and  pj = PRo^jtD, ,  ,•  Then  the  measurement  noise  covariance  matrix, 
,  for  a  sequence  of  m  pairwise  correlated  range  difference  observations  is  given  by 


P„.  = 


PlO;03  0 

^2  92^2^3 

92O2O3  O3  ^ 

0 


0 


pm  — i  ®m-i 


g 


(F7) 
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The  corresponding  linear  measurement  model  for  these  m  observations  is  given  by 

jt'=  A'Ax+V  (F8) 

where  £(v')  =  0  and  fi(v'  v  ')  =  ?„•.  To  whiten  and  decorrelate  these  observations, 
i.e.,  to  transform  the  observations  into  an  equivalent  set  of  independent  observations 
each  with  unit  variance,  do  the  following: 

Factor  P„>  into  the  form  LL'  where  L  is  a  lower  triaiuralar  matrix.  Then  transform 
the  linear  measurement  model  given  by  equation  (F8)  by  multiplying  each  side  by  L ' 


to  get 

L-’z'  =  L-'AAx+L'V  (P9) 

or  X  =  AAx+v  (FIO) 

where  E(v)  =  E(L'v' )  =  L'B(v')=0  (Fll) 

andE(w^)  =  ECLVv'a’^)  (F12) 


=  L'E(v'V^)L^ 
=  L‘P„<L^ 


=  L‘LL^L^ 


=  I 


Therefore  eouation  (FIO)  is  the  linear  measurement  model  for  the  equivalent  set  of 
uncorrelatea,  unit  variance  observations  required  for  the  Filter  measurement 
update  step. 

The  banded  tridiagonal  matrix,  ?„•,  has  a  lower  triangular  square  root,  L,  of  the 

form 


L 


Ot 

Pi  0 

Pa  03 

0 

Pm-J 


(F13) 


Then  is  also  given  by 
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P„.  =  LL^  = 


'  8i  fli  Pf  0 

t  Jt 

it  I  Pi  Og+Pi  8aPi 

0  OgPi 


Om-I  +  Pm-«  ®m-iPm-J 
Om-lPm-1  8m  +  Pm-1 


(P14) 


Equating  terms  between  the  matrices  given  in  equations  (F14)  and  (F7)  gives 


5,  =  Oi 

8jPi  =  PiOjOg  -*  Pi  = 


PlOiOg 


oa  +  pi  —  Oj 


8j  —  (oa— Pj)* 


(F16) 


it 3  Pa  ~  PaOgOa  -»  Pa  = 


PaOaOa 


,  _i  »  _  -  t  t.  ^ 

03+P2  =  Oa  -♦  03  —  lOa— pa) 


or,  in  general, 


PnOnOn+l 


Pn  = 


On+1  —  (®n+|— P«)* 


n  =  1,  2, ...»  m-i 

8i  =  Oi 


(F16) 


The  L  ‘  matrix  is  not  actually  computed  but  the  transformation  given  by  equation 
(F9)  is  done  recursively  as  follows: 

Let  c  represent  z' ,  v ,  or  any  column  of  A.  Then  c is  given  by 


Pi  fla 


pa  83 


'  Pm-I  8,, 


(F17) 
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Individual  equations  are  then  given  by 


Ci 

9,  Cl 

Cs 

s 

PiCl 

+  9iCi 

Cj 

P«Ci 

+  93  Cj 

Cm 

s 

Pm-|Cm-J 

+  8m  Cm 

Solving  those  equations  recursively  for  the  Ci*s  then  gives 

C|  =  cj/fl/ 

=  (ci— PiCil/Oj 
=  (ci— PjCjI/Oj 


F.6 


(F18) 


(F19) 
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APPENDIX  G 

HOUSEHOLDER  ORTHOGONAL  TRANSFORMATIONS 
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This  appendix  is  not  intended  to  mvide  a  detailed  explanation  of  Householder 
transformationa  and  their  propertiee.  like  reader  ia  referred  to  ChMter  IV  of  refer¬ 
ence  4  for  this  explanation.  However*  two  properties  of  theae  transKurmationa  ex¬ 
ploited  in  the  MsF/S  implementation  will  be  meeuseed  belovr .  First  it  will  be  ahown 
that  the  tranafoimation  defined  by  equationa  (lOS-110)  in  the  MEASUREMENT 
UPDATE  aubaection  of  the  FILTER  ALGORITHM  ae^on  doea  aero  out  the  below 
diagonal  elementa  in  the  first  column  of  the  arbitrary  m  Xn  matrix  R. 

For  the  first  column,  equation  (109)  gives 


Squaring  equation  ( 105)  gives 

s»  =  |Jl(l,J)|*  (G2) 

Substituting  equation  (G2)  in  (Gl)  gives 


Yj 


—  (a»-aR(J.J))=  — (s-R(J,J)) 
ai(i)  u(l) 


(G3) 


Substituting  in  equation  (G3)  for  u(l)  from  equation  (106)  gives 
a-R(J,I) 

Yi  =  -  =  -1 

RiUD-s 

Therefore  for  J=l,  equation  (110)  gives 


(G4) 


(G5) 


Two  properties  of  Householder  orthomnal  transformations  are  exploited  to 
reduce  the  number  of  computations  in  the  MSF/S  implementation.  These  are  given 
as  follows: 


1.  If  the  current  column  being  leroed  out  below  the  diagonal  has  a  zero  ele¬ 
ment,  the  corresponding  row  for  allremaining  columns  is  unchani^  by  this  trans¬ 
formation.  u(0 = 0  for  this  row  in  equation  ( 1 10)  so  that  R(i,y)  for  all  columns  /does 
not  chanj^.  This  property  allows  the  Householder  transformations  required  for  the 
measurement  up<ute  to  be  done  in  two  stages  (equations  (103)  and  (104))  and  saves 
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over  Ny  X  (Np  +  N,)  words  of  array  storage.  Since  the  Ry  array  is  always  upper  trian- 
^lar,  this  property  allows  it  to  be  storedand  operated  on  as  a  one-dimensional  array 
in  equation  (104).  The  Filter  propagation  step  computations  are  also  reduced  by 
applying  this  property.  In  equation  (111)  the  matrix  -  Ru,M  is  always  upper  trian¬ 
gular  so  that  only  the  highest  row  being  operated  on  in  the  first  Np  rows  changes  as 
each  column  is  zeroed  below  the  diagonal. 

2.  If  any  column  being  transformed  has  zero  entries  corresponding  to  all  non¬ 
zero  entries  in  the  current  column  being  zeroed  out  below  the  diagonal,  this  entire 
column  is  unchanged  by  this  transformation,  yj  =  0  for  this  column  in  equa¬ 
tion  (110)  so  that  R(iJ)  for  all  rows  i  does  not  change.  This  property  is  applied  as 
described  under  equation  (111)  in  the  PROPAGATION  subsection  of  Uie  FILTER 
ALGORITHM  section.  The  Gauss-Markov  p  parameters  were  placed  first  in  the  list 
of  stochastic  parameters  in  order  to  take  maximum  advantage  of  this  property. 


NSWCTR  87-187 


APPENDIX  H 

SRIF  AND  EQUIVALENCE  TO  KALMAN  FILTER 
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The  square  root  information  filter  (SRIF)  measurement  update  equations  (equa¬ 
tions  (103)  and  (104))  corresponding  to  the  parameter  set  partiuoning  implemented  in 
the  MSF/S  system  are  derived  in  reference  4,  section  VII.2  and  are  pven  there  by  equa¬ 
tions  (2.6)  and  (2.7).  The  SRIF  propagation  equations  (equation  (111))  are  derived  in 
reference  4,  section  Vn.3  and  are  pven  there  In  a  more  general  form  by  equation  (3.6) 
repeated  here  as 


Rp  Rpp  Rpx  ^py 


Hp—  RpjVp  0  Rpx  Rpy—RpxVy  Zp 

Rxp~Rx^p  ®  Rx  R*y~Rx^v 


Rp  Rpx  ^py  ^p 


^xp  Rx  ^ 


xy 


(HI) 


where  Rpx=  RpxVx  and  A,  =  A,  V,  (H2) 

These  equations  reduce  to  that  used  in  the  MSF/S  system  with  the  following 
simplifications: 

Zu,  =  0  (H3) 

i.e.,  no  a  priori  knowledge  of  process  noise  exists  except  for  its  covariance  matrix, 

Vy  =  0  (H4) 

and  Vj,  =  i  (H5) 

because  of  the  pseudoepoch  state  variable  definitions  (see  Appendices  A  and  D),  and 

Rxp  =  0  (H6) 


because  the  measurement  update  Householder  transformation  is  applied  even  if  no 
observations  are  present  in  a  given  mini-batch  interval. 

An  alternate  derivation  of  the  propagation  equations  can  be  obtained  by 
substituting  partitioned  ^  and  G  matrices  into  the  general  propagation  algorithm  given 
by  equation  (2.29)  in  reference  4,  section  VI.2  and  rearranging  the  result  as  follows: 

The  general  form  of  the  state  equations  is  given  by  equation  (1)  and  repeated 
here  as 


Axy+/  =  ^j^Xj  GWj 


(H7) 


The  MSF/S  state  equations  correspond  to  this  form  under  the  following  definitions: 


^XJ 


(H8) 
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(H9) 


(HIO) 


The  general  form  of  the  informaticn  array  being  transformed  in  the  propagation  step  is 
given  by  the  following,  taken  from  reference  4,  equation  (2.29),  with  rjj)  =  0: 


/  RJi)  0  0  \ 

I  -I  iL  ,  I  *  / 

\-Rj^jG  Rj<pj  Xj/ 

This  is  equivalent  to  the  data  equations  given  by 
RwU)'^j  = 

-Rj^‘jG}gj  +  Rj^j^Xj  =  Zj-Vj 

For  the  partitioning  done  in  the  MSF/S  system  and  ij  are  given  by 


Inverting  both  sides  of  equation  (H9)  gives 


Then  Rj^j  and  -Rj^C  are  given  by 

'kpM'-Rp^WpM'  Rp, 


-R^VpM‘  R, 


(HID 

(H12) 

(H13) 


(H14) 


(H15) 


(H16) 


(H17) 


H.4 


0 


0 


NSWCTR  87-187 


-(R„AI'-Jlp,VpM') 


—  Rj^jG  — 


RxVpM' 


(H18) 


Substituting  equations  (H15),  (H17K  nnd  (H18)  into  equation  (H13)  gives 
-(RpM  '-Rp,VpM  ')w^  +  (Rp^M  '-Rp,VpM  ')Ap,+,  +  Rp,^  Axj>/ +  Rp,  Ay,+,  =  ip  -  Vp^  (H19) 
KVpM  '  wj -  R,^ Vp, M  •  Apj + ,  +  R,^  Axj + /  +  R,,  Ay, + ,  =  i,^ -  (H20) 

R,Ay,.,  =i,-v,^(H21) 


Solving  the  state  equations  for  Wj  gives 

Mfj  =  —  MAp,  +  Apj+ j  (H22) 

Multiplying  each  side  of  equation  (H22)  by  Rj/)  gives 

RJi)yfj=  -Ru,(j)MAp,-l-  R  J;)Ap,+i  (H23) 

Ignoring  equation  (H21),  since  y  parameters  do  not  change  in  the  propagation  step,  and 
substituting  equation  (H23)  into  equation  (H12)  and  equation  (H22)  into  equations 
(H19)  and  (H20)  gives 

-Ru,(/)MAp,  +  Ru,0)Apj^.;  =  (H24) 

(Rp^ -Rp^Vp^Apj  -|-Rp,^Ax,>,  +Rpy^Ax,+,  =  ip,  -Vp^  (H25) 

-R,Vp^  Ap,  +Rx,  Ax,+  ,  +R,y^  Ay,>,  =  i,,  -v,^  (H26) 


This  is  a  set  of  data  equations  for  the  states  Ap,,  Apy  + 1,  Axj+;,  and  corresponds 

to  the  information  array  required  for  the  propagation  step  in  the  MSF/S  implementation 
given  by 


(H27) 


The  Householder  transformation,  Tp,  operating  on  this  matrix  eliminates  Ap/  from  the 
last  Np -I- N.  rows. 


The  SRIF  algorithm  is  mathematically  equivalent  to  the  standard  Kalmar  filter 
algorithm.  Assume  a  set  of  state  equations  given  by  equation  (1)  (repeated  above  as 
equation  (H7))  and  a  linear  measurement  model  given  by  equation  (2)  and  repeated  here 
as 


Zj  =  A^Ax^-l-y, 


(H28) 
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Then  this  equivalence  will  be  shown  for  each  filter  processing  step-measurement 
update  and  propagation. 

Measurement  Update 

The  Kalman  filter  measurement  update  equations  are  given  by 

dxj  =  S^Xj-k-Kj{Zj  —  Aj£xj)  (H29) 

Pj  =  (/  -  KjAj)Pj  =  ( Pj'+  A  ■  Aj)  ‘  (H30) 

Kj  =  PjAj(AjPjAj  +  n  =  Pj  Aj  (H31) 

where  '  indicates  a  predi9ted  quantity  and  ''  indicates  a  filter  estimated  quantity.  A 
different  expression  for  Ax^  can  be  derived  by  substituting  for  Kj  from  equation  (H31) 
into  equation  (H29)  and  using  equation  (H30)  to  give 

AXj  =  ^Xj  +  PjA](2j-Aj£Xj) 

=  (I-Pj  AjAj)Sxj  +  PjAjZj 

-  r  *  .  r  -  r  ,  (H32) 

=  Pj[(Pj-AjAj)£iXj-¥AjZj] 

~  Pj(Pj'Sxj+  aJzj) 

The  SRIF  measurement  update  equations  are  given  in  non- partitioned  form  by 


(H33) 


Premultiply  each  side  by  its  transpose  to  get 


.r «  .  r 
ZjZj+9j0j 


Since  f  s  orthogonal,  t’f =/  so  equation  (H34)  reduces  to 

/-r--  r  •  7‘  \ 

/RjRj+  AjAj  R,2j+AjzA  ^  /RjRj  RjZj 

I  J  »  \  .r«  ,  r 

\IjRj+  ZjAj  ^7^7+  \f/^y 

Equating  upper  left  terms  gives 

a]Aj  =  k]Rj 

Inverting  both  sides  gives 

RjkJ  =  (ftjftj+AjAj)' 


(H34) 


(H35) 


(H36) 


(H37) 
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Subtituting  P  =  R '  A  ^  in  equation  (H37)  then  gives 

Pj  =  (Pj+  /LjA/  (H38) 

which  is  the  same  as  equation  (M30)  above. 

Equating  upper  right  terms  in  equation  (H35)  gives 

FiiZj+ AjZj  =  (H39) 

Mulitplying  z  =  RAxby  R^  on  each  side  gives 

R^z  =  R’RAx  (H40) 


Substituting  equation  (H40)  into  equation  (H39)  for  terms  involving  "  and '  gives 

rJRjAIj  +  A]zj  =  rJRjAXj  (H41) 

Since  P  =  R  'R  ^  P '  =  R^R  so  that  equation  (H41)  becomes 

p/A*Xj  +  AjZj  =  PjAxj  (H42) 

Multiplying  each  side  by  Pj  and  regrouping  terms  then  gives 

Axj  =  Pj  (Pj^Xj  +  A^Zj)  (H43) 

which  is  the  same  as  equation  (H32)  above. 

Propagation 

The  Kalman  filter  propagation  equations  are  given  by 

Axj+i  =  ^jdxj  (H44) 

Pj^j  =  ^jPj^j+GQjC'^  (H45) 

The  SRIF  propagation  equations  are  given  by 


RJi)  0  0 

y—  Rj^j  C  Rj<Pj  ij 
Partition  f  as  follows: 


Kii)  RM  iwii) 

0  Pj+i  Ij  +  i 


(H46) 


f 


fn 

Ta, 


} 

}  N, 


(H47) 
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Using  these  definitions  it  follows  that 

ftj+i  =  or  ^22=  (H48) 

f'ilRu){j)  —  f22^j^jG  =0  or 

.  -  .1  (H49) 

^2/  —  Ru>(/)  —  Rj+jGRttiO) 

Also,  since  f  is  orthogonal  it  follows  that 

^22^21  "I"  22^22  ~  ^  (H50) 

Substituting  for  T23  &nd  72i  from  equations  (H48)  and  (H49)  into  equation  (H50)  gives 
I  =  RtoiJ)  Rw(.0  G  Rj+i  +  Rj  Rj^jRj^j  (H51) 

Pre-multiplying  by  Rj\  1  and  post-multiplying  by  RJ^  i  gives 

i  ^ + 1  =  C  r1(/)  r1(/)  G  +  Rj  Rj%  (H52) 

Substituting  P  =  R '  R  ^and  Qj  =  rU})  RLH)  then  gives 

Pj+i  =  ^jPj^]+CQjG'  (H53) 

which  is  the  same  as  equation  (H45)  above. 

Also  the  definitions  in  equations  (H46)  and  (H47)  give 

f22ij  =  ij  (H54) 

Substituting  for  ^22  from  equ  tion  (H48)  above  gives 

Rj^i^j  r‘zj  =  2j+,  (H55) 

Since  2  =  rAx,  this  becomes 

=  Rj^iH^xj+t  (H56) 

Multiplying  both  sides  by  j  then  gives 

^jdxj  =  Axj+/  (H57) 

which  is  the  same  as  equation  (H44)  above. 
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APPENDIX! 

SRIS  AND  EQUIVALENCE  TO  RTS  SMOOTHER 
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The  square  root  information  smoother  (SRIS)  is  implemented  using  two  algo¬ 
rithms  in  the  MSF/S  system.  The  first  alTOrithm,  called  the  State  Only  SmooUier. 
provides  state  estimates  only.  The  second  algorithm,  called  the  Array  Smoother,  in 
addition  provides  covariance  estimates.  Both  smoothing  algorithms  are  developed  in 
reference  4.  section  X.2  for  the  general  form  of  the  state  equations  given  by  equation 
(H7).  For  the  special  foim  of  the  state  equations  defined  by  equations  (H8),  (H9),  and 
(HIO),  the  Array  Smoother  and  State  Only  Smoother  can  be  derived  by  the  proper 
substitutions  as  follows: 


The  information  array,  corresponding  to  the  general  form  of  the  state  equa¬ 
tions,  that  is  transformed  in  the  SRiS  algorithm  is  given  in  reference  4,  section  X.2, 
equation  2.7  and  is  repeated  here  as 


ftj/)  +  RuJi)c 

Rxitj+i)C 


(ID 


The  data  equation  saved  for  smoothing  in  the  general  case  based  on  reference  4, 
section  VI.2,  equation  2.29  is  given  by 


RJjlWj  +  =  iuij  -  ^u>j  (12) 

The  corresponding  data  equation  saved  for  smoothing  in  the  MSF/S  propagation  step 
is  given  by 


^p(tj)A|%  Spp(tj)A|i,+  ,  +  ftp,(tj)A*^+i  +  ftpy(tj)Ayjv  ~  fp^  -  Vp^  (13) 

From  the  state  equations  it  follows  that 


Apj.y;  =  MA|^  +  Wj  or  Apj  M''Ap^4)  -  M  'w, 

Substituting  for  A|i^  in  equation  (13)  above  gives 

—  Rp(tj}M ’Mfj  ■4-  '  +  /lpp(tj))Apf -f }  4*  Apx(ty)A]^^j  +  Rpyitj)^ys  ~  ~^Pj 

Comparing  terms  between  equations  (15)  and  (12)  gives 
Rjj)  =  -Rp(tj)M' 

Rwxij)  =  (ftp(tj)AI'  Rppitj)  Rpxitj)  Rpy[tj)) 

Therefore 

RJi)  +  RMC  =  -Sp(tj)M'-l-ftp(tj)M'  +  ftpp(tj)  =  Rppitj) 

RM^j  =  (Rp(tj) +Rpp(tj)M+ Rpx{tj)Vp^  Rpxitj)  Rpyitj)) 


(14) 

(15) 

(16) 

(17) 

(18) 
(19) 

(110) 
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Also 


0  Rlitj^,)  Riyitj^i)  I  I  0  j  =  I 

0  0  R;{tj^,)  /  \.0  J  \ 

^R}>(tj^i)  RM^,)  Rpy{tj^,)\  /m  0  0 

0  RKtj^i)  Riyitj^,)  I  Vp^  /  0 

K  0  0  R;(tj.,)/  \0  0  I 

(Rpitj-i-  l)M  +  Rpi(tj+|)Vp,  Rpyi^j-¥l) 

0  0  R*yitj^,) 

—  I  +  I 


(Ill) 


(112) 


i*yitj  +  l) 


(113) 


Substituting  orations  (18)  through  (113)  into  equation  (II)  gives  the  sm< 
array  for  the  Nl(SF/S  implementation  as 

!  Rppltj)  Rpitj)  +  RppUj)M  +  ftpx(tj)Vp,  RpAtj)  Rpyitj)  ip^ 

I  Rpitj-^l)  Rpitj-t-  l)M  +  Rpxi^j  +  l)ypj  R^(tj+l)  Rpy(^J*  l)  l^pity'+z) 


the  smoothing 


R:(tj+/)Vp^ 


Riitj+i)  R:,(tj^,)  ziitj^,) 

0  R*y(tj^l)  x;(tj  +  |) 


(114) 


Also  R*y{tj)  =  ftj,(tjv)  and  z^(tj)  =  Zy(t^)  for  all  /  =  0, 1, N  where  Rj,(tw)  is  upper 
triangular.  These  last  Ny  rows  are  therefore  not  carried  in  the  smoothing  compu¬ 
tations.  The  matrices  ftp(t;),  Rp(tj^  i),  and  Rldj^i)  are  always  upper  triangular.  The 
Householder  transformation,  Tpp^,  operating  on  this  matrix  eliminates  from 

the  last  Np  +  N,  +  Ny  rows. 

The  State  Only  Smoother  algorithm  corresponding  to  the  general  form  of  the 
state  equations  is  given  in  reference  4,  section  X.2,  equations  (2.1)  and  (2.2)  and  is 
repeated  he.  las 


Ax^  =  RsZn 

Wj  =  [Sj/)]'(?J/)-^u.x(/)Ax;+j) 

Ax;  =  $)'(Ax;+i-cw*) 


i  =  N-1,  N~2,...,0 


(115) 


(116) 
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Solving  equation  (13)  for  Ap/by  setting  =  0  gives 

Apr  =  (117) 

This  is  the  same  as  equation  (114)  in  the  STATE  ONLY  SMOOTHER  subsection  of  the 
SMOOTHER  Al^GORITHMS  section.  To  derive  equation  (115)  substitute  the 

partitioned  definitions  of  Ax,  end  G  into  the  second  part  of  equation  (116)  to  get 


M 


=  l-VpM' 


(118) 


Extracting  out  the  portion  of  equation  (118)  involving  Axj* gives 

Ax;  =  -Vp^Ai'(Ap;+;-w;)+Ax;+i  (Ii9) 

Prom  the  state  equations 

Apr+;  =  AlApr+w/or  yfj-  Af^Vi-MAp,*  (120) 

Substituting  w/into  equation  (119)  then  gives 
Ax;  =  -VpM'(MAp;)+Ax;+, 

=  Ax;+/-Vp^Ap; 

This  is  the  same  as  equation  (115).  Therefore  equations  (115),  (117),  and  (121)  give 
the  State  Only  Smoother  algorithm  in  the  MSF^  i 


(121) 


i  system. 


The  SRIS  al^rithm  is  mathematically  equivalent  to  the  standard  Rauch- 
Tung-Striebel  (RTS)  smoothing  algorithm.  Assume  a  set  of  state  equations  of  the 
general  form  given  by  equation  (H7).  Then  this  equivalence  can  be  shown  by  proper 
manipulation  of  the  smoother  information  arrays  as  follows: 


The  RTS  smoother  equations  are  given  by 
Ax;  =  dxj  +  Cj{  Ax;+ ,  -  Axj + 1 ) 

P]  =  Pj  +  CjiP]^t-Pj*i)C/ 
where  Cj  =  Pj^jPj'+j 

The  SRIS  equations  are  given  by 

y  Rj\,G  RJ^,%  t%i  j 


RM 


(122) 

(123) 

(124) 


(125) 
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Using  these  dv^finitions  it  follows  that 

0  =  +ThR;->-tG 

Rj  =  ThRwxmj  +ThRUi% 

Post-multiplying  equation  (128)  by  G  gives 
=  ThRu,xU)G  +n2Rj^iG 

Subtracting  equation  (127)  from  equation  (129)  gives 
Rji^'jG  =  -ThRwij)  or  Th  -  -R;4>jc  CO) 


(126) 

(127) 

(128) 

(129) 


(130) 


Post-multiplying  equation  (128)  by  4>>  and  substituting  for  from  equation  (130) 
gives 


Rj^j  =  -RJ^jGRMRmW  +  ThRUl 


(131) 


Solving  this  for  Tls  gives 

Th  =  RjI^^I+GrUDRwxW)  Rj'ii 

In  reference  4,  section  X.A  it  is  shown  from  tine  propagation  equations  that 
^j(l+GRUj}Ru,x()))  =  Pj^jP%i  =  Cj 
Therefore 

rp*  _  D*  ^ 

1 22  —  RjGjHj  +  j 

Since  r*is  orthogonal  it  follows  that 

fW9*  _L  W 

^21^21  'T  I  22  122  —  i 


Substituting  for  Th  and  Th  from  equations  (130)  and  (134)  gives 


/  =  RJi)  RJi)  G  (^jPTj  +R]Cj  fCj^ifCj^iCjR) 

Pre-multiplying  by  R^'and  post-multiplying  by  fC/sind  replacing 
Kd)  rLU)  by  Qj ,  Rj+iRj+t  by  P;+/,  and  ff/  R^^  by  p; gives 
~  GQjG  +  CjPj  + ;  Cj 


(132) 

(133) 

(134) 

(135) 

(136) 

(137) 
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In  reference  4,  section  X.A  it  is  also  shown  that 
i>jbQjG%^=  Pj-Pj  ^]pUi^jPj 

Using  the  second  part  of  equation  (133)  in  equation  (138)  gives 

=  Pj-Cj^jPj 

Substituting  equation  (139)  into  equation  (137)  gives 
p;  =  Pj-CjPj^,cf+CjP;^,cJ 
Rearranging  terms  then  gives 

p;  =  Pj+Cj(p;+i-Pj^i)cl 
which  is  the  same  as  equation  (123)  above. 

Also  the  definitions  in  equations  (125)  and  (126)  give 

=  T'2lZj.j)  +  T22Zj  +  i 

Substituting  for  Th  and  T22  from  equations  (130)  and  (134)  give 
z;  =  -R*^jGRUimi)+RjCjRf^,z;^, 

Pre-multi plying  by  R^and  substituting  Ax*  =  R*' z*  gives 
Ax;  =  -^jG  R‘ji)zJi)  +  Cj^x;^i 
In  reference  4,  section  X.A  it  is  also  shown  that 

—  ^^G  Ru;(/)Zu)(7)  —  Aj^— C^Aj^+i 

Substituting  this  into  equation  (144)  then  gives 
Ax;  =  iSxj  +  C,(Ax;+;-AX;  +  ;) 
which  is  the  same  as  equation  (122)  above. 


(138) 

(139) 

(140) 

(141) 

(142) 

(143) 

(144) 

(145) 

(146) 
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APPENDIX  J 

INVERSION  OP  UPPERTRIANGULAR  MATRICES 
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Let  If  be  a  nonainmlar  upper  triangular  matrix  of  dimension  n  X  n  to  be 
inverted  In  the  MSF/S  implementation,  the  lower  triangular  matrix  L = U'^  is 
actually  computed  to  take  maximum  advantage  of  existing  array  storage  space.  Let 
Uij  t  U  and  tijtL,  then  L  is  computed  as  follows: 

Define 

=  l/ui,t  (Jl) 

For  each  subsequent  row  j  from  j  =  2, ....  n  define 
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